NS_P-delta_(v0.8) -- Sage

%auto

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


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

def bind_1st (f, a) :
    return lambda *x : f (a, *x)
    
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 map2 (f, *l) :
    return map (bind_1st (map, f), *l)  

     
def imax (l) :
    ix = max (xrange (len(l)), key = lambda i : l[i])
    return ix, l[ix]

def imax_ (l, f) :
    ix = max (xrange (len(l)), key = lambda i : f (l[i]))
    return ix, l[ix]

def iamax (l) :
    return imax_ (l, abs)


def extract_from (ll, els = None, ix = 0) :
    if els == None :
        els = range (0, len (ll))
    lf = []
    for ei in els :
        lf.append (ll[ei][ix]);
    return lf


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

def print_indexed (l, indent = '') :
    for i in xrange (len(l)) :
        print '{0}{1:>3}: {2}'.format (indent, 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 (tol = None) :
    def f_ (x) :
        t = sqrt ((parent(x)(1)).ulp()) if tol == None else tol
        return x if (abs(x) >= t) else parent(x)(0)
    return f_

def NDC (sig_digits_10, tolerance = None) :
    return lambda x : N (chop(tolerance) (x), digits = sig_digits_10)


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

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


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

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


##########################
## 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) :
    return map (apply_to_ith (lambda x : -x, -1), nodes)
    
def coords_left_to_right_t (nodes) :
    return map (apply_to_ith_t (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)   



##########################    
## 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, n_dof_per_node = 0) :

    if n_dof_per_node == 0 :  
        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)
    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()
    


############################    
## 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)
    


##############################    
## 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 frame2D_element_geometric_stiffness_matrix_ (hh, ll) :
    k = zero_matrix (parent (ll), 6, 6)
    k[1,1] = 6/5
    k[1,2] = k[2,1] = -ll/10  
    k[1,4] = k[4,1] = -k[1,1]
    k[1,5] = k[5,1] = k[1,2]
    k[2,2] = 2*ll^2/15
    k[2,4] = k[4,2] = -k[1,2]
    k[2,5] = k[5,2] = -ll^2/30
    k[4,4] = k[1,1]
    k[4,5] = k[5,4] = -k[1,2]
    k[5,5] = k[2,2]
    return hh/ll*k

def frame2D_element_geometric_stiffness_matrix (ch, nc) :
    (hh) = ch
    ll = element2D_length (nc)
    return frame2D_element_geometric_stiffness_matrix_ (hh, 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 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



##############################    
## 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 distributed_force2D_geometric_ (qz, hh, ee, ii, ll) :
    f = vector (parent(ll), 6)
    f[2] = -hh*qz*ll^4 / (720*ee*ii)
    f[5] = -f[2]
    return f
    
def distributed_force2D_geometric (qq, ll) :
    if len(qq) == 3 :
        qz, hh, (ee, aa, ii) = qq
    elif len(qq) == 4 :
        qz, hh, ee, ii = qq
    else :
        print "wrong element load"
        return
    return distributed_force2D_geometric_ (qz, hh, ee, ii, ll)
    
 
def fixed_end_force_vector (f, ff, enc) :
    ll = element2D_length (enc)
    return f (ff, ll)


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

def structure_stiffness_matrix (fk, fr, ndof, tdof, nodes, elems, echs, chars) :
    field = parent (ith_node_coords (0, nodes)[0])
    sk = 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
    

    
##############################    
## 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, els, nds, 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, els)
        enc = element_node_coords (el, nds)
        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 (ndof, tdof, nodal_loads = [], elem_loads = [], els = [], nds = [], fr = None) :
    if nodal_loads != [] :
        sfn = structure_force_vector_nodes (ndof, tdof, nodal_loads)
    if elem_loads != [] :
        sfe = structure_force_vector_elems (ndof, tdof, elem_loads, els, nds, fr)
    if nodal_loads != [] and elem_loads != [] :
        sf = sfn + sfe
    elif nodal_loads != [] :
        sf = sfn
    elif elem_loads != [] :
        sf = sfe
    return sf
    
    
    
##############################    
## 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.rows()
    
def nodal_displacements (sk, sf, tdof, uf = []) :
    d = sk \ 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()
    nd = len (uu[0])
    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, nds, els, elchs, chars, elem_loads = []) :
    efs = []
    for j in xrange (len (els)) :
        el = ith_element (j, els)
        enc = element_node_coords (el, nds)
        ec = ith_element_characteristics (j, elchs, 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, els, nds)
            efl = fixed_end_force_vector (f, ldv, enc)
            efs[ei] += efl

    return efs
    
    
    
##########################
##  Main function (DiM):
##########################
    
    
def DiM (nodes, elems, supports, chars, el_ch_ids, fk, fr,
         nodal_loads = [], elem_loads = [], field = RR) :

    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)

    n, t = nodal_DsOF (len (nds2), supports)
    sk = structure_stiffness_matrix (fk, fr, n, t, nds2, els, el_ch_ids, chrs2)
    sf = structure_force_vector (n, t, ndlds2, ellds2, els, nds2, fr)
    d = nodal_displacements (sk, sf, t)
    ef = all_end_forces (fk, fr, d, nds2, els, el_ch_ids, chrs2, ellds2)

    return (d, 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, cf_scale, df_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, df_scale, color = color))
        elif f.func_name == concentrated_force2D.func_name :
            pls.append (plot_concentrated (v, ec, cf_scale, color = color, width = width))
    return sum (pls)
    
    
def plot_frame (nodes, elems, supports = [], hinges = [], releases = [],
                nodal_loads = [], elem_loads = [], cf_scale = 0.01, df_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, cf_scale)
    if elem_loads != [] :
        pt += plot_elem_loads (elem_loads, nodes, elems, cf_scale, df_scale)
    return pt

%auto

# Kratice:

f2k = frame2D_element_stiffness_matrix
f2r = frame2D_element_global_to_local
f2kg = frame2D_element_geometric_stiffness_matrix
df2 = distributed_force2D
df2g = distributed_force2D_geometric

nds_xy = [(0., 0.), (-5., 4.), (0., 4.), (5., 4.)]
print 'nodes:'
print_list (nds_xy, '   ')
print 

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

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

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

chars = [(3e7, 100000, 0.3^4/12), (3e7, 100000, 2*0.3^4/12)]
print 'characteristics (large axial stiffness):'
print_list (chars, '   ')
print

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

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

ellds = [(2, df2, (0., 50.))]
print 'element loads:'
print_list (ellds, '   ')

nodes:
   (0.000000000000000, 0.000000000000000)
   (-5.00000000000000, 4.00000000000000)
   (0.000000000000000, 4.00000000000000)
   (5.00000000000000, 4.00000000000000)

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

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

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

characteristics (large axial stiffness):
   (3.00000000000000e7, 100000, 0.000675000000000000)
   (3.00000000000000e7, 100000, 0.00135000000000000)

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

nodal loads:
   (2, (0.000000000000000, 1000.00000000000, 0.000000000000000))
   (3, (250.000000000000, 0.000000000000000, 0.000000000000000))

element loads:
   (2, <function distributed_force2D at 0x7f8e9a3e8320>,
(0.000000000000000, 50.0000000000000))

nodes:
   (0.000000000000000, 0.000000000000000)
   (-5.00000000000000, 4.00000000000000)
   (0.000000000000000, 4.00000000000000)
   (5.00000000000000, 4.00000000000000)

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

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

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

characteristics (large axial stiffness):
   (3.00000000000000e7, 100000, 0.000675000000000000)
   (3.00000000000000e7, 100000, 0.00135000000000000)

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

nodal loads:
   (2, (0.000000000000000, 1000.00000000000, 0.000000000000000))
   (3, (250.000000000000, 0.000000000000000, 0.000000000000000))

element loads:
   (2, <function distributed_force2D at 0x7f8e9a3e8320>, (0.000000000000000, 50.0000000000000))

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds, cf_scale = 0.0025, df_scale = 0.01)

nds = coords_left_to_right (nds_xy)
print_list (nds)

[0.000000000000000, -0.000000000000000]
[-5.00000000000000, -4.00000000000000]
[0.000000000000000, -4.00000000000000]
[5.00000000000000, -4.00000000000000]

[0.000000000000000, -0.000000000000000]
[-5.00000000000000, -4.00000000000000]
[0.000000000000000, -4.00000000000000]
[5.00000000000000, -4.00000000000000]

ndof, tbl_dof = nodal_DsOF (4, supps)
ndof; tbl_dof

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

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

# opća metoda pomaka 
K0 = structure_stiffness_matrix (f2k, f2r, ndof, tbl_dof, nds, els, echs, chars0)
show (K0.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 32400.0 & 0.000000 & 9720.00 & 16200.0 & 0.000000 & 0.000000 \\ 0.000000 & 2.16380 \times 10^{6} & 0.000000 & 7593.75 & -1.08000 \times 10^{6} & 0.000000 \\ 9720.00 & 0.000000 & 682776. & 0.000000 & 0.000000 & -9720.00 \\ 16200.0 & 7593.75 & 0.000000 & 85050.0 & 0.000000 & 16200.0 \\ 0.000000 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 & 1.08000 \times 10^{6} & 0.000000 \\ 0.000000 & 0.000000 & -9720.00 & 16200.0 & 0.000000 & 32400.0 \end{array}\right)$

# "inženjerska" metoda pomaka
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}{rrrrrr} 32400.0 & 0.000000 & 9720.00 & 16200.0 & 0.000000 & 0.000000 \\ 0.000000 & 1.20000 \times 10^{12} & 0.000000 & 7593.75 & -6.00000 \times 10^{11} & 0.000000 \\ 9720.00 & 0.000000 & 7.50000 \times 10^{11} & 0.000000 & 0.000000 & -9720.00 \\ 16200.0 & 7593.75 & 0.000000 & 85050.0 & 0.000000 & 16200.0 \\ 0.000000 & -6.00000 \times 10^{11} & 0.000000 & 0.000000 & 6.00000 \times 10^{11} & 0.000000 \\ 0.000000 & 0.000000 & -9720.00 & 16200.0 & 0.000000 & 32400.0 \end{array}\right)$

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

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.000000,\,0.000000,\,1125.00,\,-104.167,\,250.000,\,104.167\right)$

u0 = nodal_displacements (K0, q, tbl_dof)
print_indexed (map2 (ND(6), u0))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.000635908]
  2: [0.000246762, 0.00170804, -0.00229664]
  3: [0.000478244, 0.000000, 0.00487575]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.000635908]
  2: [0.000246762, 0.00170804, -0.00229664]
  3: [0.000478244, 0.000000, 0.00487575]

u = nodal_displacements (K, q, tbl_dof)
print_indexed (map2 (ND(6), u))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.00113471]
  2: [4.45389e-10, 1.54167e-9, -0.00226943]
  3: [8.62056e-10, 0.000000, 0.00434973]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.00113471]
  2: [4.45389e-10, 1.54167e-9, -0.00226943]
  3: [8.62056e-10, 0.000000, 0.00434973]

print_indexed (map2 (NDC(6), u))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.00113471]
  2: [0.000000, 0.000000, -0.00226943]
  3: [0.000000, 0.000000, 0.00434973]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.00113471]
  2: [0.000000, 0.000000, -0.00226943]
  3: [0.000000, 0.000000, 0.00434973]

ef0 = all_end_forces (f2k, f2r, u0, nds, els, echs, chars0, ellds)
print_indexed (map2 (ND(6), ef0))

  0: [1152.93, 16.5032, -21.3797, -1152.93, -16.5032, -44.6332]
  1: [-266.503, 9.50146, -7.10543e-15, 266.503, -9.50146, -47.5073]
  2: [-250.000, -143.428, 92.1405, 250.000, -106.572, 0.000000]

  0: [1152.93, 16.5032, -21.3797, -1152.93, -16.5032, -44.6332]
  1: [-266.503, 9.50146, -7.10543e-15, 266.503, -9.50146, -47.5073]
  2: [-250.000, -143.428, 92.1405, 250.000, -106.572, 0.000000]

print_indexed (map2 (NDC(6), ef0))

  0: [1152.93, 16.5032, -21.3797, -1152.93, -16.5032, -44.6332]
  1: [-266.503, 9.50146, 0.000000, 266.503, -9.50146, -47.5073]
  2: [-250.000, -143.428, 92.1405, 250.000, -106.572, 0.000000]

  0: [1152.93, 16.5032, -21.3797, -1152.93, -16.5032, -44.6332]
  1: [-266.503, 9.50146, 0.000000, 266.503, -9.50146, -47.5073]
  2: [-250.000, -143.428, 92.1405, 250.000, -106.572, 0.000000]

ef = all_end_forces (f2k, f2r, u, nds, els, echs, chars, ellds)
print_indexed (map2 (NDC(6), ef))

  0: [1156.25, 17.2335, -22.9779, -1156.25, -17.2335, -45.9559]
  1: [-267.233, 11.0294, 0.000000, 267.233, -11.0294, -55.1471]
  2: [-250.000, -145.221, 101.103, 250.000, -104.779, 0.000000]

  0: [1156.25, 17.2335, -22.9779, -1156.25, -17.2335, -45.9559]
  1: [-267.233, 11.0294, 0.000000, 267.233, -11.0294, -55.1471]
  2: [-250.000, -145.221, 101.103, 250.000, -104.779, 0.000000]

print_indexed (map2 (lambda x, y : x - y, u0, u))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, -0.000498804490057595]
  2: [0.000246761781700479, 0.00170804225516144, -0.0000272163729816708]
  3: [0.000478242846515294, 0.000000000000000, 0.000526020863039265]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, -0.000498804490057595]
  2: [0.000246761781700479, 0.00170804225516144, -0.0000272163729816708]
  3: [0.000478242846515294, 0.000000000000000, 0.000526020863039265]

print_indexed (map2 (lambda x, y : x - y, ef0, ef))

  0: [-3.32043414403392, -0.730249307564442, 1.59828150334859, 3.32043414403392,
0.730249307564442, 1.32271572690917]
  1: [0.730249307564463, -1.52794549932600, 0.000000000000000,
-0.730249307564463, 1.52794549932600, 7.63972749662997]
  2: [0.000000000000000, 1.79248864470784, -8.96244322353918, 0.000000000000000,
-1.79248864470783, -2.84217094304040e-14]

  0: [-3.32043414403392, -0.730249307564442, 1.59828150334859, 3.32043414403392, 0.730249307564442, 1.32271572690917]
  1: [0.730249307564463, -1.52794549932600, 0.000000000000000, -0.730249307564463, 1.52794549932600, 7.63972749662997]
  2: [0.000000000000000, 1.79248864470784, -8.96244322353918, 0.000000000000000, -1.79248864470783, -2.84217094304040e-14]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef0, ef)))

  0: [-3.32043414403392, -0.730249307564442, 1.59828150334859, 3.32043414403392,
0.730249307564442, 1.32271572690917]
  1: [0.730249307564463, -1.52794549932600, 0.000000000000000,
-0.730249307564463, 1.52794549932600, 7.63972749662997]
  2: [0.000000000000000, 1.79248864470784, -8.96244322353918, 0.000000000000000,
-1.79248864470783, 0.000000000000000]

  0: [-3.32043414403392, -0.730249307564442, 1.59828150334859, 3.32043414403392, 0.730249307564442, 1.32271572690917]
  1: [0.730249307564463, -1.52794549932600, 0.000000000000000, -0.730249307564463, 1.52794549932600, 7.63972749662997]
  2: [0.000000000000000, 1.79248864470784, -8.96244322353918, 0.000000000000000, -1.79248864470783, 0.000000000000000]

print_indexed (map2 (lambda x, y : abs (x - y)/abs (x) * 100, ef0, ef))

  0: [0.287999740053513, 4.42489380822793, 7.47571169721131, 0.287999740053513,
4.42489380822793, 2.96352667231420]
  1: [0.274011453956427, 16.0811569418188, 0.000000000000000, 0.274011453956427,
16.0811569418188, 16.0811569418188]
  2: [0.000000000000000, 1.24974720541903, 9.72693252973207, 0.000000000000000,
1.68195238013220, +infinity]

  0: [0.287999740053513, 4.42489380822793, 7.47571169721131, 0.287999740053513, 4.42489380822793, 2.96352667231420]
  1: [0.274011453956427, 16.0811569418188, 0.000000000000000, 0.274011453956427, 16.0811569418188, 16.0811569418188]
  2: [0.000000000000000, 1.24974720541903, 9.72693252973207, 0.000000000000000, 1.68195238013220, +infinity]

print_indexed (map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef0, ef))

  0: [0.287999740053513, 4.42489380822793, 7.47571169721131, 0.287999740053513,
4.42489380822793, 2.96352667231420]
  1: [0.274011453956427, 16.0811569418188, NaN, 0.274011453956427,
16.0811569418188, 16.0811569418188]
  2: [0.000000000000000, 1.24974720541903, 9.72693252973207, 0.000000000000000,
1.68195238013220, NaN]

  0: [0.287999740053513, 4.42489380822793, 7.47571169721131, 0.287999740053513, 4.42489380822793, 2.96352667231420]
  1: [0.274011453956427, 16.0811569418188, NaN, 0.274011453956427, 16.0811569418188, 16.0811569418188]
  2: [0.000000000000000, 1.24974720541903, 9.72693252973207, 0.000000000000000, 1.68195238013220, NaN]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef0, ef)))

  0: [0.287999740053513, 4.42489380822793, 7.47571169721131, 0.287999740053513,
4.42489380822793, 2.96352667231420]
  1: [0.274011453956427, 16.0811569418188, 0.000000000000000, 0.274011453956427,
16.0811569418188, 16.0811569418188]
  2: [0.000000000000000, 1.24974720541903, 9.72693252973207, 0.000000000000000,
1.68195238013220, 0.000000000000000]

  0: [0.287999740053513, 4.42489380822793, 7.47571169721131, 0.287999740053513, 4.42489380822793, 2.96352667231420]
  1: [0.274011453956427, 16.0811569418188, 0.000000000000000, 0.274011453956427, 16.0811569418188, 16.0811569418188]
  2: [0.000000000000000, 1.24974720541903, 9.72693252973207, 0.000000000000000, 1.68195238013220, 0.000000000000000]

# opća metoda pomaka 
Hs0 = extract_from (ef0, ix = 3)
print Hs0
eHs = range (len (Hs0))
print eHs

[-1152.92956285897, 266.503205256736, 250.000000000000]
[0, 1, 2]

[-1152.92956285897, 266.503205256736, 250.000000000000]
[0, 1, 2]

Kg0 = structure_stiffness_matrix (f2kg, f2r, ndof, tbl_dof, nds, els, eHs, Hs0)
show (Kg0.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 177.669 & 0.000000 & 26.6503 & -44.4172 & 0.000000 & 0.000000 \\ 0.000000 & -345.879 & 0.000000 & -115.293 & 0.000000 & 0.000000 \\ 26.6503 & 0.000000 & 123.961 & 1.65032 & 0.000000 & -25.0000 \\ -44.4172 & -115.293 & 1.65032 & -270.560 & 0.000000 & -41.6667 \\ 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & -25.0000 & -41.6667 & 0.000000 & 166.667 \end{array}\right)$

# "inženjerska" metoda pomaka
Hs = extract_from (ef, ix = 3)
print Hs

[-1156.24999700300, 267.233454564300, 250.000000000000]

[-1156.24999700300, 267.233454564300, 250.000000000000]

Kg = structure_stiffness_matrix (f2kg, f2r, ndof, tbl_dof, nds, els, eHs, Hs)
show (Kg.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 178.156 & 0.000000 & 26.7233 & -44.5389 & 0.000000 & 0.000000 \\ 0.000000 & -346.875 & 0.000000 & -115.625 & 0.000000 & 0.000000 \\ 26.7233 & 0.000000 & 124.136 & 1.72335 & 0.000000 & -25.0000 \\ -44.5389 & -115.625 & 1.72335 & -271.844 & 0.000000 & -41.6667 \\ 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & -25.0000 & -41.6667 & 0.000000 & 166.667 \end{array}\right)$

K20 = K0 + Kg0
show (K20.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 32577.7 & 0.000000 & 9746.65 & 16155.6 & 0.000000 & 0.000000 \\ 0.000000 & 2.16345 \times 10^{6} & 0.000000 & 7478.46 & -1.08000 \times 10^{6} & 0.000000 \\ 9746.65 & 0.000000 & 682900. & 1.65032 & 0.000000 & -9745.00 \\ 16155.6 & 7478.46 & 1.65032 & 84779.4 & 0.000000 & 16158.3 \\ 0.000000 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 & 1.08000 \times 10^{6} & 0.000000 \\ 0.000000 & 0.000000 & -9745.00 & 16158.3 & 0.000000 & 32566.7 \end{array}\right)$

K2 = K + Kg
show (K2.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 32578.2 & 0.000000 & 9746.72 & 16155.5 & 0.000000 & 0.000000 \\ 0.000000 & 1.20000 \times 10^{12} & 0.000000 & 7478.12 & -6.00000 \times 10^{11} & 0.000000 \\ 9746.72 & 0.000000 & 7.50000 \times 10^{11} & 1.72335 & 0.000000 & -9745.00 \\ 16155.5 & 7478.12 & 1.72335 & 84778.2 & 0.000000 & 16158.3 \\ 0.000000 & -6.00000 \times 10^{11} & 0.000000 & 0.000000 & 6.00000 \times 10^{11} & 0.000000 \\ 0.000000 & 0.000000 & -9745.00 & 16158.3 & 0.000000 & 32566.7 \end{array}\right)$

elldsg0 = [(2, df2g, (50., Hs0[2], chars0[1][0], chars0[1][2]))]
print_list (elldsg0)

(2, <function distributed_force2D_geometric at 0x7f8e9a3e8410>,
(50.0000000000000, 250.000000000000, 3.00000000000000e7, 0.00135000000000000))

(2, <function distributed_force2D_geometric at 0x7f8e9a3e8410>, (50.0000000000000, 250.000000000000, 3.00000000000000e7, 0.00135000000000000))

elldsg = [(2, df2g, (50., Hs[2], chars[1][0], chars[1][2]))]
print_list (elldsg)

(2, <function distributed_force2D_geometric at 0x7f8e9a3e8410>,
(50.0000000000000, 250.000000000000, 3.00000000000000e7, 0.00135000000000000))

(2, <function distributed_force2D_geometric at 0x7f8e9a3e8410>, (50.0000000000000, 250.000000000000, 3.00000000000000e7, 0.00135000000000000))

qg = structure_force_vector (ndof, tbl_dof, [], elldsg, els, nds, f2r)
show (qg.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.000000,\,0.000000,\,0.000000,\,0.267918,\,0.000000,\,-0.267918\right)$

q2 = q + qg
show (q2.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.000000,\,0.000000,\,1125.00,\,-103.899,\,250.000,\,103.899\right)$

u20 = nodal_displacements (K20, q2, tbl_dof)
print_indexed (map2 (ND(6), u20))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.000623767]
  2: [0.000246537, 0.00170751, -0.00228796]
  3: [0.000478018, 0.000000, 0.00483648]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.000623767]
  2: [0.000246537, 0.00170751, -0.00228796]
  3: [0.000478018, 0.000000, 0.00483648]

u2 = nodal_displacements (K2, q2, tbl_dof)
print_indexed (map2 (NDC(6), u2))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.00112127]
  2: [0.000000, 0.000000, -0.00226110]
  3: [0.000000, 0.000000, 0.00431221]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.000000, 0.000000, 0.00112127]
  2: [0.000000, 0.000000, -0.00226110]
  3: [0.000000, 0.000000, 0.00431221]

ef2e0 = all_end_forces (f2k, f2r, u20, nds, els, echs, chars0, ellds)
ef2g0 = all_end_forces (f2kg, f2r, u20, nds, els, eHs, Hs0, elldsg)
ef20 = map (lambda x, y : x + y, ef2e0, ef2g0)
print_indexed (map2 (NDC(6), ef20))

  0: [1152.57, 16.2596, -21.6736, -1152.57, -16.2596, -43.0806]
  1: [-266.260, 9.47232, 0.000000, 266.260, -9.47232, -47.8167]
  2: [-250.000, -143.094, 90.8973, 250.000, -106.906, 0.000000]

  0: [1152.57, 16.2596, -21.6736, -1152.57, -16.2596, -43.0806]
  1: [-266.260, 9.47232, 0.000000, 266.260, -9.47232, -47.8167]
  2: [-250.000, -143.094, 90.8973, 250.000, -106.906, 0.000000]

ef2e = all_end_forces (f2k, f2r, u2, nds, els, echs, chars, ellds)
ef2g = all_end_forces (f2kg, f2r, u2, nds, els, eHs, Hs, elldsg)
ef2 = map (lambda x, y : x + y, ef2e, ef2g)
print_indexed (map2 (NDC(6), ef2))

  0: [1156.10, 16.9088, -23.2422, -1156.10, -16.9088, -44.3929]
  1: [-266.909, 11.1095, 0.000000, 266.909, -11.1095, -55.5476]
  2: [-250.000, -144.988, 99.9405, 250.000, -105.012, 0.000000]

  0: [1156.10, 16.9088, -23.2422, -1156.10, -16.9088, -44.3929]
  1: [-266.909, 11.1095, 0.000000, 266.909, -11.1095, -55.5476]
  2: [-250.000, -144.988, 99.9405, 250.000, -105.012, 0.000000]

print_indexed (map2 (lambda x, y : x - y, u20, u2))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, -0.000497506765300211]
  2: [0.000246536241752720, 0.00170750424058654, -0.0000268647107679520]
  3: [0.000478017306567535, 0.000000000000000, 0.000524269737236683]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, -0.000497506765300211]
  2: [0.000246536241752720, 0.00170750424058654, -0.0000268647107679520]
  3: [0.000478017306567535, 0.000000000000000, 0.000524269737236683]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef20, ef2)))

  0: [-3.53122235092815, -0.649138249288313, 1.56857679093620, 3.53122235092815,
0.649138249288313, 1.31221512617282]
  1: [0.649138249288228, -1.63720759202999, 0.000000000000000,
-0.649138249288228, 1.63720759202999, 7.73098260816984]
  2: [0.000000000000000, 1.89401475889790, -9.04319773434271, 0.000000000000000,
-1.89401475889790, 0.000000000000000]

  0: [-3.53122235092815, -0.649138249288313, 1.56857679093620, 3.53122235092815, 0.649138249288313, 1.31221512617282]
  1: [0.649138249288228, -1.63720759202999, 0.000000000000000, -0.649138249288228, 1.63720759202999, 7.73098260816984]
  2: [0.000000000000000, 1.89401475889790, -9.04319773434271, 0.000000000000000, -1.89401475889790, 0.000000000000000]

print_indexed (map2 (lambda x, y : abs (x - y)/abs (x) * 100, ef20, ef2))

  0: [0.306379080814160, 3.99233308193736, 7.23726714376592, 0.306379080814160,
3.99233308193736, 3.04595084004958]
  1: [0.243798982948019, 17.2841255626754, 325.454545454545, 0.243798982948019,
17.2841255626754, 16.1679709879168]
  2: [6.82121026329696e-14, 1.32361500475753, 9.94880858161537,
6.82121026329696e-14, 1.77166504835466, 72.4137931034483]

  0: [0.306379080814160, 3.99233308193736, 7.23726714376592, 0.306379080814160, 3.99233308193736, 3.04595084004958]
  1: [0.243798982948019, 17.2841255626754, 325.454545454545, 0.243798982948019, 17.2841255626754, 16.1679709879168]
  2: [6.82121026329696e-14, 1.32361500475753, 9.94880858161537, 6.82121026329696e-14, 1.77166504835466, 72.4137931034483]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef20, ef2)))

  0: [0.306379080814160, 3.99233308193736, 7.23726714376592, 0.306379080814160,
3.99233308193736, 3.04595084004958]
  1: [0.243798982948019, 17.2841255626754, 0.000000000000000, 0.243798982948019,
17.2841255626754, 16.1679709879168]
  2: [0.000000000000000, 1.32361500475753, 9.94880858161537, 0.000000000000000,
1.77166504835466, 0.000000000000000]

  0: [0.306379080814160, 3.99233308193736, 7.23726714376592, 0.306379080814160, 3.99233308193736, 3.04595084004958]
  1: [0.243798982948019, 17.2841255626754, 0.000000000000000, 0.243798982948019, 17.2841255626754, 16.1679709879168]
  2: [0.000000000000000, 1.32361500475753, 9.94880858161537, 0.000000000000000, 1.77166504835466, 0.000000000000000]

print_indexed (map2 (lambda x, y : x - y, u0, u20))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000121409644914861]
  2: [2.25540488916934e-7, 5.38014778058349e-7, -8.68162767910831e-6]
  3: [2.25540488916934e-7, 0.000000000000000, 0.0000392751330329528]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000121409644914861]
  2: [2.25540488916934e-7, 5.38014778058349e-7, -8.68162767910831e-6]
  3: [2.25540488916934e-7, 0.000000000000000, 0.0000392751330329528]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, u0, u20)))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000121409644914861]
  2: [2.25540488916934e-7, 5.38014778058349e-7, -8.68162767910831e-6]
  3: [2.25540488916934e-7, 0.000000000000000, 0.0000392751330329528]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000121409644914861]
  2: [2.25540488916934e-7, 5.38014778058349e-7, -8.68162767910831e-6]
  3: [2.25540488916934e-7, 0.000000000000000, 0.0000392751330329528]

print_indexed (map2 (lambda x, y : x - y, u, u2))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000134386892488704]
  2: [5.41157977177455e-13, 2.03162357727038e-13, -8.32996546538955e-6]
  3: [5.41157977177352e-13, 0.000000000000000, 0.0000375240072303705]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000134386892488704]
  2: [5.41157977177455e-13, 2.03162357727038e-13, -8.32996546538955e-6]
  3: [5.41157977177352e-13, 0.000000000000000, 0.0000375240072303705]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, u, u2)))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000134386892488704]
  2: [0.000000000000000, 0.000000000000000, -8.32996546538955e-6]
  3: [0.000000000000000, 0.000000000000000, 0.0000375240072303705]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 0.000000000000000, 0.0000134386892488704]
  2: [0.000000000000000, 0.000000000000000, -8.32996546538955e-6]
  3: [0.000000000000000, 0.000000000000000, 0.0000375240072303705]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef0, ef20)))

  0: [0.363159975189546, 0.243583728030263, 0.293949557506853,
-0.363159975189546, -0.243583728030263, -1.55252390393908]
  1: [-0.243583728030274, 0.0291455290406990, 0.000000000000000,
0.243583728030274, -0.0291455290406990, 0.309328118707228]
  2: [0.000000000000000, -0.334014446148899, 1.24319578523185,
0.000000000000000, 0.334014446148885, 0.000000000000000]

  0: [0.363159975189546, 0.243583728030263, 0.293949557506853, -0.363159975189546, -0.243583728030263, -1.55252390393908]
  1: [-0.243583728030274, 0.0291455290406990, 0.000000000000000, 0.243583728030274, -0.0291455290406990, 0.309328118707228]
  2: [0.000000000000000, -0.334014446148899, 1.24319578523185, 0.000000000000000, 0.334014446148885, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef0, ef20)))

  0: [0.0314988865658891, 1.47597829779670, 1.37490306985351,
0.0314988865658891, 1.47597829779670, 3.47840878060776]
  1: [0.0913999243632426, 0.306747738621939, 0.000000000000000,
0.0913999243632426, 0.306747738621939, 0.651116682583609]
  2: [0.000000000000000, 0.232879366838169, 1.34923940075151, 0.000000000000000,
0.313416988362694, 0.000000000000000]

  0: [0.0314988865658891, 1.47597829779670, 1.37490306985351, 0.0314988865658891, 1.47597829779670, 3.47840878060776]
  1: [0.0913999243632426, 0.306747738621939, 0.000000000000000, 0.0913999243632426, 0.306747738621939, 0.651116682583609]
  2: [0.000000000000000, 0.232879366838169, 1.34923940075151, 0.000000000000000, 0.313416988362694, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef, ef2)))

  0: [0.152371768295325, 0.324694786306392, 0.264244845094460,
-0.152371768295325, -0.324694786306392, -1.56302450467544]
  1: [-0.324694786306509, -0.0801165636632923, 0.000000000000000,
0.324694786306509, 0.0801165636632923, 0.400583230247094]
  2: [0.000000000000000, -0.232488331958848, 1.16244127442832,
0.000000000000000, 0.232488331958820, 0.000000000000000]

  0: [0.152371768295325, 0.324694786306392, 0.264244845094460, -0.152371768295325, -0.324694786306392, -1.56302450467544]
  1: [-0.324694786306509, -0.0801165636632923, 0.000000000000000, 0.324694786306509, 0.0801165636632923, 0.400583230247094]
  2: [0.000000000000000, -0.232488331958848, 1.16244127442832, 0.000000000000000, 0.232488331958820, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef, ef2)))

  0: [0.0131780989137533, 1.88409575744037, 1.14999371022852,
0.0131780989137533, 1.88409575744037, 3.40114149886368]
  1: [0.121502297246389, 0.726390260180758, 0.000000000000000,
0.121502297246389, 0.726390260180758, 0.726391007148389]
  2: [0.000000000000000, 0.160093233036189, 1.14976019567017, 0.000000000000000,
0.221883597348537, 0.000000000000000]

  0: [0.0131780989137533, 1.88409575744037, 1.14999371022852, 0.0131780989137533, 1.88409575744037, 3.40114149886368]
  1: [0.121502297246389, 0.726390260180758, 0.000000000000000, 0.121502297246389, 0.726390260180758, 0.726391007148389]
  2: [0.000000000000000, 0.160093233036189, 1.14976019567017, 0.000000000000000, 0.221883597348537, 0.000000000000000]

nds_xy = [(0., 0.), (-5., 4.), (0., 4.), (5., 4.)]
print 'nodes:'
print_list (nds_xy, '   ')
print 

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

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

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

chars = [(3e7, 100000, 0.3^4/12), (3e7, 100000, 2*0.3^4/12)]
print 'characteristics (large axial stiffness):'
print_list (chars, '   ')
print

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

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

ellds = [(2, df2, (0., 50.))]
print 'element loads:'
print_list (ellds, '   ')

nodes:
   (0.000000000000000, 0.000000000000000)
   (-5.00000000000000, 4.00000000000000)
   (0.000000000000000, 4.00000000000000)
   (5.00000000000000, 4.00000000000000)

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

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

characteristics:
   (3.00000000000000e7, 100000, 0.000675000000000000)
   (3.00000000000000e7, 100000, 0.00135000000000000)

characteristics (large axial stiffness):
   (3.00000000000000e7, 100000, 0.000675000000000000)
   (3.00000000000000e7, 100000, 0.00135000000000000)

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

nodal loads:
   (2, (0.000000000000000, 1000.00000000000, 0.000000000000000))
   (3, (250.000000000000, 0.000000000000000, 0.000000000000000))

element loads:
   (2, <function distributed_force2D at 0x7f8e9a3e8320>,
(0.000000000000000, 50.0000000000000))

nodes:
   (0.000000000000000, 0.000000000000000)
   (-5.00000000000000, 4.00000000000000)
   (0.000000000000000, 4.00000000000000)
   (5.00000000000000, 4.00000000000000)

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

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

characteristics:
   (3.00000000000000e7, 100000, 0.000675000000000000)
   (3.00000000000000e7, 100000, 0.00135000000000000)

characteristics (large axial stiffness):
   (3.00000000000000e7, 100000, 0.000675000000000000)
   (3.00000000000000e7, 100000, 0.00135000000000000)

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

nodal loads:
   (2, (0.000000000000000, 1000.00000000000, 0.000000000000000))
   (3, (250.000000000000, 0.000000000000000, 0.000000000000000))

element loads:
   (2, <function distributed_force2D at 0x7f8e9a3e8320>, (0.000000000000000, 50.0000000000000))

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds, cf_scale = 0.0025, df_scale = 0.01)

nds = coords_left_to_right (nds_xy)
print_list (nds)

[0.000000000000000, -0.000000000000000]
[-5.00000000000000, -4.00000000000000]
[0.000000000000000, -4.00000000000000]
[5.00000000000000, -4.00000000000000]

[0.000000000000000, -0.000000000000000]
[-5.00000000000000, -4.00000000000000]
[0.000000000000000, -4.00000000000000]
[5.00000000000000, -4.00000000000000]

ndof, tbl_dof = nodal_DsOF (4, supps)
ndof; tbl_dof

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

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

# "inženjerska" metoda pomaka
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} 6.00000 \times 10^{11} & 0.000000 & -6.00000 \times 10^{11} & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 32400.0 & 0.000000 & 9720.00 & 16200.0 & 0.000000 & 0.000000 \\ -6.00000 \times 10^{11} & 0.000000 & 1.20000 \times 10^{12} & 0.000000 & 7593.75 & -6.00000 \times 10^{11} & 0.000000 \\ 0.000000 & 9720.00 & 0.000000 & 7.50000 \times 10^{11} & 0.000000 & 0.000000 & -9720.00 \\ 0.000000 & 16200.0 & 7593.75 & 0.000000 & 85050.0 & 0.000000 & 16200.0 \\ 0.000000 & 0.000000 & -6.00000 \times 10^{11} & 0.000000 & 0.000000 & 6.00000 \times 10^{11} & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -9720.00 & 16200.0 & 0.000000 & 32400.0 \end{array}\right)$

# opća metoda pomaka
K0 = structure_stiffness_matrix (f2k, f2r, ndof, tbl_dof, nds, els, echs, chars0)
show (K0.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 1.08000 \times 10^{6} & 0.000000 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 32400.0 & 0.000000 & 9720.00 & 16200.0 & 0.000000 & 0.000000 \\ -1.08000 \times 10^{6} & 0.000000 & 2.16380 \times 10^{6} & 0.000000 & 7593.75 & -1.08000 \times 10^{6} & 0.000000 \\ 0.000000 & 9720.00 & 0.000000 & 682776. & 0.000000 & 0.000000 & -9720.00 \\ 0.000000 & 16200.0 & 7593.75 & 0.000000 & 85050.0 & 0.000000 & 16200.0 \\ 0.000000 & 0.000000 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 & 1.08000 \times 10^{6} & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -9720.00 & 16200.0 & 0.000000 & 32400.0 \end{array}\right)$

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

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.000000,\,0.000000,\,0.000000,\,1125.00,\,-104.167,\,250.000,\,104.167\right)$

u = nodal_displacements (K, q, tbl_dof)
print_indexed (map2 (NDC(6), u))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0903020, 0.000000, 0.00611460]
  2: [0.0903020, 0.000000, -0.0122292]
  3: [0.0903020, 0.000000, 0.00932963]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0903020, 0.000000, 0.00611460]
  2: [0.0903020, 0.000000, -0.0122292]
  3: [0.0903020, 0.000000, 0.00932963]

u0 = nodal_displacements (K0, q, tbl_dof)
print_indexed (map2 (NDC(6), u0))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0903020, 0.000000, 0.00560219]
  2: [0.0903020, 0.00170804, -0.0122292]
  3: [0.0905335, 0.000000, 0.00984204]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0903020, 0.000000, 0.00560219]
  2: [0.0903020, 0.00170804, -0.0122292]
  3: [0.0905335, 0.000000, 0.00984204]

ef = all_end_forces (f2k, f2r, u, nds, els, echs, chars, ellds)
print_indexed (map2 (NDC(6), ef))

  0: [1156.25, -250.000, 561.910, -1156.25, 250.000, 438.090]
  1: [-7.62939e-6, 59.4340, 0.000000, 7.62939e-6, -59.4340, -297.170]
  2: [-250.000, -96.8160, -140.920, 250.000, -153.184, 0.000000]

  0: [1156.25, -250.000, 561.910, -1156.25, 250.000, 438.090]
  1: [-7.62939e-6, 59.4340, 0.000000, 7.62939e-6, -59.4340, -297.170]
  2: [-250.000, -96.8160, -140.920, 250.000, -153.184, 0.000000]

u[1][0]; u[2][0]; u[1][0] - u[2][0]

0.0903020376259434
0.0903020376259434
-1.38777878078145e-17

0.0903020376259434
0.0903020376259434
-1.38777878078145e-17

chars = [(3e7, 10000, 0.3^4/12), (3e7, 10000, 2*0.3^4/12)]
print 'characteristics (large axial stiffness):'
print_list (chars, '   ')
print

characteristics (large axial stiffness):
   (3.00000000000000e7, 10000, 0.000675000000000000)
   (3.00000000000000e7, 10000, 0.00135000000000000)

characteristics (large axial stiffness):
   (3.00000000000000e7, 10000, 0.000675000000000000)
   (3.00000000000000e7, 10000, 0.00135000000000000)

# "inženjerska" metoda pomaka (ponovno)
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} 6.00000 \times 10^{10} & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 32400.0 & 0.000000 & 9720.00 & 16200.0 & 0.000000 & 0.000000 \\ -6.00000 \times 10^{10} & 0.000000 & 1.20000 \times 10^{11} & 0.000000 & 7593.75 & -6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 9720.00 & 0.000000 & 7.50000 \times 10^{10} & 0.000000 & 0.000000 & -9720.00 \\ 0.000000 & 16200.0 & 7593.75 & 0.000000 & 85050.0 & 0.000000 & 16200.0 \\ 0.000000 & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -9720.00 & 16200.0 & 0.000000 & 32400.0 \end{array}\right)$

u = nodal_displacements (K, q, tbl_dof)
print_indexed (map2 (NDC(6), u))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0903020, 0.000000, 0.00611460]
  2: [0.0903020, 1.54167e-8, -0.0122292]
  3: [0.0903020, 0.000000, 0.00932963]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0903020, 0.000000, 0.00611460]
  2: [0.0903020, 1.54167e-8, -0.0122292]
  3: [0.0903020, 0.000000, 0.00932963]

u[1][0]; u[2][0]; u[1][0] - u[2][0]

0.0903020419692863
0.0903020419692863
0.000000000000000

0.0903020419692863
0.0903020419692863
0.000000000000000

ef = all_end_forces (f2k, f2r, u, nds, els, echs, chars, ellds)
print_indexed (map2 (NDC(6), ef))

  0: [1156.25, -250.000, 561.910, -1156.25, 250.000, 438.090]
  1: [0.000000, 59.4339, 0.000000, 0.000000, -59.4339, -297.170]
  2: [-250.000, -96.8160, -140.920, 250.000, -153.184, 0.000000]

  0: [1156.25, -250.000, 561.910, -1156.25, 250.000, 438.090]
  1: [0.000000, 59.4339, 0.000000, 0.000000, -59.4339, -297.170]
  2: [-250.000, -96.8160, -140.920, 250.000, -153.184, 0.000000]

ef0 = all_end_forces (f2k, f2r, u0, nds, els, echs, chars0, ellds)
print_indexed (map2 (NDC(6), ef0))

  0: [1152.93, -250.000, 561.910, -1152.93, 250.000, 438.090]
  1: [0.000000, 57.7737, 0.000000, 0.000000, -57.7737, -288.869]
  2: [-250.000, -95.1558, -149.221, 250.000, -154.844, 0.000000]

  0: [1152.93, -250.000, 561.910, -1152.93, 250.000, 438.090]
  1: [0.000000, 57.7737, 0.000000, 0.000000, -57.7737, -288.869]
  2: [-250.000, -95.1558, -149.221, 250.000, -154.844, 0.000000]

print_indexed (map2 (lambda x, y : x - y, u0, u))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [1.14729267752800e-10, 0.000000000000000, -0.000512408507721568]
  2: [1.14729267752800e-10, 0.00170802838016183, -1.26539646816415e-11]
  3: [0.000231477429544089, 0.000000000000000, 0.000512408520375530]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [1.14729267752800e-10, 0.000000000000000, -0.000512408507721568]
  2: [1.14729267752800e-10, 0.00170802838016183, -1.26539646816415e-11]
  3: [0.000231477429544089, 0.000000000000000, 0.000512408520375530]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef0, ef)))

  0: [0.287997400535099, 1.35808659251775e-7, 1.32246022489517e-7,
0.287997400535099, 1.35808659251775e-7, 1.40378271306100e-7]
  1: [0.000000000000000, 2.87362981499488, 0.000000000000000, 0.000000000000000,
2.87362981499488, 2.87362981499486]
  2: [0.000000000000000, 1.74472109176110, 5.56290573432979, 0.000000000000000,
1.07217697046569, 0.000000000000000]

  0: [0.287997400535099, 1.35808659251775e-7, 1.32246022489517e-7, 0.287997400535099, 1.35808659251775e-7, 1.40378271306100e-7]
  1: [0.000000000000000, 2.87362981499488, 0.000000000000000, 0.000000000000000, 2.87362981499488, 2.87362981499486]
  2: [0.000000000000000, 1.74472109176110, 5.56290573432979, 0.000000000000000, 1.07217697046569, 0.000000000000000]

# "inženjerska" metoda pomaka (ponovno)
Hs = extract_from (ef, ix = 3)
print Hs
eHs = range (len (Hs))
print eHs

[-1156.24997003000, 0.000000000000000, 250.000000000000]
[0, 1, 2]

[-1156.24997003000, 0.000000000000000, 250.000000000000]
[0, 1, 2]

# "opća" metoda pomaka (ponovno)
Hs0 = extract_from (ef0, ix = 3)
print Hs0

[-1152.92956285897, 0.000000000000000, 250.000000000000]

[-1152.92956285897, 0.000000000000000, 250.000000000000]

Kg = structure_stiffness_matrix (f2kg, f2r, ndof, tbl_dof, nds, els, eHs, Hs)
show (Kg.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 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 \\ 0.000000 & 0.000000 & -346.875 & 0.000000 & -115.625 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & 60.0000 & -25.0000 & 0.000000 & -25.0000 \\ 0.000000 & 0.000000 & -115.625 & -25.0000 & -450.000 & 0.000000 & -41.6667 \\ 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -25.0000 & -41.6667 & 0.000000 & 166.667 \end{array}\right)$

Kg0 = structure_stiffness_matrix (f2kg, f2r, ndof, tbl_dof, nds, els, eHs, Hs0)
show (Kg0.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 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 \\ 0.000000 & 0.000000 & -345.879 & 0.000000 & -115.293 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & 60.0000 & -25.0000 & 0.000000 & -25.0000 \\ 0.000000 & 0.000000 & -115.293 & -25.0000 & -448.229 & 0.000000 & -41.6667 \\ 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -25.0000 & -41.6667 & 0.000000 & 166.667 \end{array}\right)$

K2 = K + Kg
show (K2.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 6.00000 \times 10^{10} & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 32400.0 & 0.000000 & 9720.00 & 16200.0 & 0.000000 & 0.000000 \\ -6.00000 \times 10^{10} & 0.000000 & 1.20000 \times 10^{11} & 0.000000 & 7478.12 & -6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 9720.00 & 0.000000 & 7.50000 \times 10^{10} & -25.0000 & 0.000000 & -9745.00 \\ 0.000000 & 16200.0 & 7478.12 & -25.0000 & 84600.0 & 0.000000 & 16158.3 \\ 0.000000 & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -9745.00 & 16158.3 & 0.000000 & 32566.7 \end{array}\right)$

K20 = K0 + Kg0
show (K20.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 1.08000 \times 10^{6} & 0.000000 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 32400.0 & 0.000000 & 9720.00 & 16200.0 & 0.000000 & 0.000000 \\ -1.08000 \times 10^{6} & 0.000000 & 2.16345 \times 10^{6} & 0.000000 & 7478.46 & -1.08000 \times 10^{6} & 0.000000 \\ 0.000000 & 9720.00 & 0.000000 & 682836. & -25.0000 & 0.000000 & -9745.00 \\ 0.000000 & 16200.0 & 7478.46 & -25.0000 & 84601.8 & 0.000000 & 16158.3 \\ 0.000000 & 0.000000 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 & 1.08000 \times 10^{6} & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -9745.00 & 16158.3 & 0.000000 & 32566.7 \end{array}\right)$

elldsg = [(2, df2g, (50., Hs[2], chars[1][0], chars[1][2]))]
print_list (elldsg)

(2, <function distributed_force2D_geometric at 0x7f8e9a3e8410>,
(50.0000000000000, 250.000000000000, 3.00000000000000e7, 0.00135000000000000))

(2, <function distributed_force2D_geometric at 0x7f8e9a3e8410>, (50.0000000000000, 250.000000000000, 3.00000000000000e7, 0.00135000000000000))

qg = structure_force_vector (ndof, tbl_dof, [], elldsg, els, nds, f2r)
qg.apply_map (ND(6))

(0.000000, 0.000000, 0.000000, 0.000000, 0.267918, 0.000000, -0.267918)

(0.000000, 0.000000, 0.000000, 0.000000, 0.267918, 0.000000, -0.267918)

q2 = q + qg
q2.apply_map (ND(6))

(0.000000, 0.000000, 0.000000, 1125.00, -103.899, 250.000, 103.899)

(0.000000, 0.000000, 0.000000, 1125.00, -103.899, 250.000, 103.899)

u2 = nodal_displacements (K2, q2, tbl_dof)
print_indexed (map2 (NDC(6), u2))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.101380, 0.000000, 0.00667013]
  2: [0.101380, 1.54057e-8, -0.0133403]
  3: [0.101380, 0.000000, 0.00980928]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.101380, 0.000000, 0.00667013]
  2: [0.101380, 1.54057e-8, -0.0133403]
  3: [0.101380, 0.000000, 0.00980928]

u20 = nodal_displacements (K20, q2, tbl_dof)
print_indexed (map2 (NDC(6), u20))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.101340, 0.000000, 0.00615541]
  2: [0.101340, 0.00170667, -0.0133348]
  3: [0.101571, 0.000000, 0.0103173]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.101340, 0.000000, 0.00615541]
  2: [0.101340, 0.00170667, -0.0133348]
  3: [0.101571, 0.000000, 0.0103173]

ef2e = all_end_forces (f2k, f2r, u2, nds, els, echs, chars, ellds)
ef2g = all_end_forces (f2kg, f2r, u2, nds, els, eHs, Hs, elldsg)
ef2 = map (lambda x, y : x + y, ef2e, ef2g)
print_indexed (map2 (NDC(6), ef2))

  0: [1155.42, -250.000, 621.004, -1155.42, 250.000, 496.217]
  1: [0.000000, 64.8337, 0.000000, 0.000000, -64.8337, -324.169]
  2: [-250.000, -90.5904, -172.048, 250.000, -159.410, 0.000000]

  0: [1155.42, -250.000, 621.004, -1155.42, 250.000, 496.217]
  1: [0.000000, 64.8337, 0.000000, 0.000000, -64.8337, -324.169]
  2: [-250.000, -90.5904, -172.048, 250.000, -159.410, 0.000000]

ef2e0 = all_end_forces (f2k, f2r, u20, nds, els, echs, chars0, ellds)
ef2g0 = all_end_forces (f2kg, f2r, u20, nds, els, eHs, Hs0, elldsg)
ef20 = map (lambda x, y : x + y, ef2e0, ef2g0)
print_indexed (map2 (NDC(6), ef20))

  0: [1152.00, -250.000, 620.802, -1152.00, 250.000, 496.036]
  1: [0.000000, 63.1484, 0.000000, 0.000000, -63.1484, -315.742]
  2: [-250.000, -88.8559, -180.294, 250.000, -161.144, 0.000000]

  0: [1152.00, -250.000, 620.802, -1152.00, 250.000, 496.036]
  1: [0.000000, 63.1484, 0.000000, 0.000000, -63.1484, -315.742]
  2: [-250.000, -88.8559, -180.294, 250.000, -161.144, 0.000000]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, u20, u2)))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.0000397671151294854, 0.000000000000000, -0.000514717063708150]
  2: [-0.0000397671151294854, 0.00170665757800130, 5.43958061551865e-6]
  3: [0.000191710199685322, 0.000000000000000, 0.000507988235645147]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.0000397671151294854, 0.000000000000000, -0.000514717063708150]
  2: [-0.0000397671151294854, 0.00170665757800130, 5.43958061551865e-6]
  3: [0.000191710199685322, 0.000000000000000, 0.000507988235645147]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef20, ef2)))

  0: [-3.41983757609751, -6.36538231901795e-7, -0.201916482656770,
3.41983757609751, 6.36538231901795e-7, -0.180551790098889]
  1: [0.000000000000000, -1.68530752760870, 0.000000000000000,
0.000000000000000, 1.68530752760870, 8.42653763804344]
  2: [0.000000000000000, 1.73453004848895, -8.24598584794447, 0.000000000000000,
-1.73453004848895, 0.000000000000000]

  0: [-3.41983757609751, -6.36538231901795e-7, -0.201916482656770, 3.41983757609751, 6.36538231901795e-7, -0.180551790098889]
  1: [0.000000000000000, -1.68530752760870, 0.000000000000000, 0.000000000000000, 1.68530752760870, 8.42653763804344]
  2: [0.000000000000000, 1.73453004848895, -8.24598584794447, 0.000000000000000, -1.73453004848895, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef20, ef2)))

  0: [0.296859801917635, 2.54615292760695e-7, 0.0325251105781384,
0.296859801917635, 2.54615292760695e-7, 0.0363989243613772]
  1: [0.000000000000000, 2.66880459463875, 0.000000000000000, 0.000000000000000,
2.66880459463875, 2.66880459463874]
  2: [0.000000000000000, 1.95207166292034, 4.57363207596644, 0.000000000000000,
1.07638418622539, 0.000000000000000]

  0: [0.296859801917635, 2.54615292760695e-7, 0.0325251105781384, 0.296859801917635, 2.54615292760695e-7, 0.0363989243613772]
  1: [0.000000000000000, 2.66880459463875, 0.000000000000000, 0.000000000000000, 2.66880459463875, 2.66880459463874]
  2: [0.000000000000000, 1.95207166292034, 4.57363207596644, 0.000000000000000, 1.07638418622539, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, u, u2)))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.0110777343372989, 0.000000000000000, -0.000555531387740787]
  2: [-0.0110777343372989, 0.000000000000000, 0.00111106276887463]
  3: [-0.0110777343372989, 0.000000000000000, -0.000479647358620917]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.0110777343372989, 0.000000000000000, -0.000555531387740787]
  2: [-0.0110777343372989, 0.000000000000000, 0.00111106276887463]
  3: [-0.0110777343372989, 0.000000000000000, -0.000479647358620917]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef, ef2)))

  0: [0.825868486108448, -2.97001918170281e-7, -59.0933728746616,
-0.825868486108448, 2.97001918170281e-7, -58.1269892534680]
  1: [0.000000000000000, -5.39976506743394, 0.000000000000000,
0.000000000000000, 5.39976506743394, 26.9988253371697]
  2: [0.000000000000000, -6.22563355354242, 31.1281639162984, 0.000000000000000,
6.22563355354242, 0.000000000000000]

  0: [0.825868486108448, -2.97001918170281e-7, -59.0933728746616, -0.825868486108448, 2.97001918170281e-7, -58.1269892534680]
  1: [0.000000000000000, -5.39976506743394, 0.000000000000000, 0.000000000000000, 5.39976506743394, 26.9988253371697]
  2: [0.000000000000000, -6.22563355354242, 31.1281639162984, 0.000000000000000, 6.22563355354242, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef, ef2)))

  0: [0.0714264655148072, 1.18800767429450e-7, 10.5165121225564,
0.0714264655148072, 1.18800767429450e-7, 13.2682872021895]
  1: [0.000000000000000, 9.08532130241883, 0.000000000000000, 0.000000000000000,
9.08532130241883, 9.08532130241884]
  2: [0.000000000000000, 6.43037523407594, 22.0892627815029, 0.000000000000000,
4.06415453354882, 0.000000000000000]

  0: [0.0714264655148072, 1.18800767429450e-7, 10.5165121225564, 0.0714264655148072, 1.18800767429450e-7, 13.2682872021895]
  1: [0.000000000000000, 9.08532130241883, 0.000000000000000, 0.000000000000000, 9.08532130241883, 9.08532130241884]
  2: [0.000000000000000, 6.43037523407594, 22.0892627815029, 0.000000000000000, 4.06415453354882, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, u0, u20)))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.0110379671074401, 0.000000000000000, -0.000553222831754204]
  2: [-0.0110379671074401, 1.37081317210593e-6, 0.00110562317560514]
  3: [-0.0110379671074401, 0.000000000000000, -0.000475227073890534]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.0110379671074401, 0.000000000000000, -0.000553222831754204]
  2: [-0.0110379671074401, 1.37081317210593e-6, 0.00110562317560514]
  3: [-0.0110379671074401, 0.000000000000000, -0.000475227073890534]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef0, ef20)))

  0: [0.925298891171451, 0.000000000000000, -58.8914556489007,
-0.925298891171451, 0.000000000000000, -57.9464368483865]
  1: [0.000000000000000, -5.37466106384428, 0.000000000000000,
0.000000000000000, 5.37466106384428, 26.8733053192214]
  2: [0.000000000000000, -6.29995995501579, 31.0731315291650, 0.000000000000000,
6.29995995501579, 0.000000000000000]

  0: [0.925298891171451, 0.000000000000000, -58.8914556489007, -0.925298891171451, 0.000000000000000, -57.9464368483865]
  1: [0.000000000000000, -5.37466106384428, 0.000000000000000, 0.000000000000000, 5.37466106384428, 26.8733053192214]
  2: [0.000000000000000, -6.29995995501579, 31.0731315291650, 0.000000000000000, 6.29995995501579, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef0, ef20)))

  0: [0.0802563244953967, 0.000000000000000, 10.4805780462258,
0.0802563244953967, 0.000000000000000, 13.2270736063069]
  1: [0.000000000000000, 9.30294753330371, 0.000000000000000, 0.000000000000000,
9.30294753330371, 9.30294753330371]
  2: [0.000000000000000, 6.62067754791726, 20.8235781047593, 0.000000000000000,
4.06858037612818, 0.000000000000000]

  0: [0.0802563244953967, 0.000000000000000, 10.4805780462258, 0.0802563244953967, 0.000000000000000, 13.2270736063069]
  1: [0.000000000000000, 9.30294753330371, 0.000000000000000, 0.000000000000000, 9.30294753330371, 9.30294753330371]
  2: [0.000000000000000, 6.62067754791726, 20.8235781047593, 0.000000000000000, 4.06858037612818, 0.000000000000000]

nds_xy = [(0., 0.), (-5., 4.), (0., 4.), (5., 4.)]
print 'nodes:'
print_list (nds_xy, '   ')
print 

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

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

chars = [(3e7, 10000, 0.000675)]
print 'characteristics:'
print_list (chars, '   ')
print

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

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

ellds = [(1, df2, (0., 40.)), (2, df2, (0., 50.))]
print 'element loads:'
print_list (ellds, '   ')

nodes:
   (0.000000000000000, 0.000000000000000)
   (-5.00000000000000, 4.00000000000000)
   (0.000000000000000, 4.00000000000000)
   (5.00000000000000, 4.00000000000000)

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

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

characteristics:
   (3.00000000000000e7, 10000, 0.000675000000000000)

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

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

element loads:
   (1, <function distributed_force2D at 0x7f8e9a3e8320>,
(0.000000000000000, 40.0000000000000))
   (2, <function distributed_force2D at 0x7f8e9a3e8320>,
(0.000000000000000, 50.0000000000000))

nodes:
   (0.000000000000000, 0.000000000000000)
   (-5.00000000000000, 4.00000000000000)
   (0.000000000000000, 4.00000000000000)
   (5.00000000000000, 4.00000000000000)

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

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

characteristics:
   (3.00000000000000e7, 10000, 0.000675000000000000)

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

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

element loads:
   (1, <function distributed_force2D at 0x7f8e9a3e8320>, (0.000000000000000, 40.0000000000000))
   (2, <function distributed_force2D at 0x7f8e9a3e8320>, (0.000000000000000, 50.0000000000000))

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds, cf_scale = 0.003, df_scale = 0.01)

nds = coords_left_to_right (nds_xy)
print_list (nds)

[0.000000000000000, -0.000000000000000]
[-5.00000000000000, -4.00000000000000]
[0.000000000000000, -4.00000000000000]
[5.00000000000000, -4.00000000000000]

[0.000000000000000, -0.000000000000000]
[-5.00000000000000, -4.00000000000000]
[0.000000000000000, -4.00000000000000]
[5.00000000000000, -4.00000000000000]

ndof, tbl_dof = nodal_DsOF (4, supps)
ndof; tbl_dof

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

7
[-1 -1 -1]
[ 0 -1  1]
[ 2  3  4]
[ 5 -1  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} 6.00000 \times 10^{10} & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 16200.0 & 0.000000 & 4860.00 & 8100.00 & 0.000000 & 0.000000 \\ -6.00000 \times 10^{10} & 0.000000 & 1.20000 \times 10^{11} & 0.000000 & 7593.75 & -6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 4860.00 & 0.000000 & 7.50000 \times 10^{10} & 0.000000 & 0.000000 & -4860.00 \\ 0.000000 & 8100.00 & 7593.75 & 0.000000 & 52650.0 & 0.000000 & 8100.00 \\ 0.000000 & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -4860.00 & 8100.00 & 0.000000 & 16200.0 \end{array}\right)$

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

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(100.000,\,-83.3333,\,0.000000,\,1225.00,\,-20.8333,\,0.000000,\,104.167\right)$

u = nodal_displacements (K, q, tbl_dof)
print_indexed (map2 (NDC(6), u))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0420888, 0.000000, -0.00120619]
  2: [0.0420888, 1.70833e-8, -0.00787569]
  3: [0.0420888, 0.000000, 0.0103679]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0420888, 0.000000, -0.00120619]
  2: [0.0420888, 1.70833e-8, -0.00787569]
  3: [0.0420888, 0.000000, 0.0103679]

ef = all_end_forces (f2k, f2r, u, nds, els, echs, chars, ellds)
print_indexed (map2 (NDC(6), ef))

  0: [1281.25, -100.000, 239.871, -1281.25, 100.000, 160.129]
  1: [100.000, -55.8621, 0.000000, -100.000, -144.138, -220.690]
  2: [0.000000, -137.112, 60.5603, 0.000000, -112.888, 0.000000]

  0: [1281.25, -100.000, 239.871, -1281.25, 100.000, 160.129]
  1: [100.000, -55.8621, 0.000000, -100.000, -144.138, -220.690]
  2: [0.000000, -137.112, 60.5603, 0.000000, -112.888, 0.000000]

Hs = extract_from (ef, ix = 3)
print Hs
eHs = range (len (Hs))
print eHs

[-1281.24998339500, -100.000000476837, 0.000000000000000]
[0, 1, 2]

[-1281.24998339500, -100.000000476837, 0.000000000000000]
[0, 1, 2]

Kg = structure_stiffness_matrix (f2kg, f2r, ndof, tbl_dof, nds, els, eHs, Hs)
show (Kg.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & -66.6667 & 0.000000 & -10.0000 & 16.6667 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & -384.375 & 0.000000 & -128.125 & 0.000000 & 0.000000 \\ 0.000000 & -10.0000 & 0.000000 & -24.0000 & -10.0000 & 0.000000 & 0.000000 \\ 0.000000 & 16.6667 & -128.125 & -10.0000 & -750.000 & 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 & 0.000000 & 0.000000 \end{array}\right)$

K2 = K + Kg
show (K2.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 6.00000 \times 10^{10} & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ 0.000000 & 16133.3 & 0.000000 & 4850.00 & 8116.67 & 0.000000 & 0.000000 \\ -6.00000 \times 10^{10} & 0.000000 & 1.20000 \times 10^{11} & 0.000000 & 7465.62 & -6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 4850.00 & 0.000000 & 7.50000 \times 10^{10} & -10.0000 & 0.000000 & -4860.00 \\ 0.000000 & 8116.67 & 7465.62 & -10.0000 & 51900.0 & 0.000000 & 8100.00 \\ 0.000000 & 0.000000 & -6.00000 \times 10^{10} & 0.000000 & 0.000000 & 6.00000 \times 10^{10} & 0.000000 \\ 0.000000 & 0.000000 & 0.000000 & -4860.00 & 8100.00 & 0.000000 & 16200.0 \end{array}\right)$

elldsg = [(1, df2g, (40., Hs[1], chars[0])), (2, df2g, (50., Hs[2], chars[0]))]
elldsg

[(1,
  <function distributed_force2D_geometric at 0x7f8e9a3e8410>,
  (40.0000000000000,
   -100.000000476837,
   (3.00000000000000e7, 10000, 0.000675000000000000))),
 (2,
  <function distributed_force2D_geometric at 0x7f8e9a3e8410>,
  (50.0000000000000,
   0.000000000000000,
   (3.00000000000000e7, 10000, 0.000675000000000000)))]

[(1,
  <function distributed_force2D_geometric at 0x7f8e9a3e8410>,
  (40.0000000000000,
   -100.000000476837,
   (3.00000000000000e7, 10000, 0.000675000000000000))),
 (2,
  <function distributed_force2D_geometric at 0x7f8e9a3e8410>,
  (50.0000000000000,
   0.000000000000000,
   (3.00000000000000e7, 10000, 0.000675000000000000)))]

qg = structure_force_vector (ndof, tbl_dof, [], elldsg, els, nds, f2r)
show (qg.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.000000,\,-0.171468,\,0.000000,\,0.000000,\,0.171468,\,0.000000,\,0.000000\right)$

q2 = q + qg
show (q2.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(100.000,\,-83.5048,\,0.000000,\,1225.00,\,-20.6619,\,0.000000,\,104.167\right)$

u2 = nodal_displacements (K2, q2, tbl_dof)
print_indexed (map2 (NDC(6), u2))

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0492012, 0.000000, -0.000600291]
  2: [0.0492012, 1.70823e-8, -0.00909489]
  3: [0.0492012, 0.000000, 0.0109775]

  0: [0.000000, 0.000000, 0.000000]
  1: [0.0492012, 0.000000, -0.000600291]
  2: [0.0492012, 1.70823e-8, -0.00909489]
  3: [0.0492012, 0.000000, 0.0109775]

ef2e = all_end_forces (f2k, f2r, u2, nds, els, echs, chars, ellds)
ef2g = all_end_forces (f2kg, f2r, u2, nds, els, eHs, Hs, elldsg)
ef2 = map (lambda x, y : x + y, ef2e, ef2g)
print_indexed (map2 (NDC(6), ef2))

  0: [1281.17, -100.000, 273.678, -1281.17, 100.000, 189.361]
  1: [100.000, -52.9784, 0.000000, -100.000, -147.022, -235.108]
  2: [0.000000, -134.149, 45.7471, 0.000000, -115.851, 0.000000]

  0: [1281.17, -100.000, 273.678, -1281.17, 100.000, 189.361]
  1: [100.000, -52.9784, 0.000000, -100.000, -147.022, -235.108]
  2: [0.000000, -134.149, 45.7471, 0.000000, -115.851, 0.000000]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, u, u2)))

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.00711234953928723, 0.000000000000000, -0.000605900829033263]
  2: [-0.00711234953928724, 0.000000000000000, 0.00121919740180058]
  3: [-0.00711234953928724, 0.000000000000000, -0.000609598700584369]

  0: [0.000000000000000, 0.000000000000000, 0.000000000000000]
  1: [-0.00711234953928723, 0.000000000000000, -0.000605900829033263]
  2: [-0.00711234953928724, 0.000000000000000, 0.00121919740180058]
  3: [-0.00711234953928724, 0.000000000000000, -0.000609598700584369]

print_indexed (map2 (chop(), map2 (lambda x, y : x - y, ef, ef2)))

  0: [0.0789797357153930, 6.13838153640245e-7, -33.8074190956349,
-0.0789797357153930, -6.13838153640245e-7, -29.2315964718309]
  1: [4.76837158203125e-7, -2.88366995014816, 0.000000000000000,
-4.76837158203125e-7, 2.88366995014817, 14.4183480425128]
  2: [0.000000000000000, -2.96264968586362, 14.8132484293181, 0.000000000000000,
2.96264968586362, 0.000000000000000]

  0: [0.0789797357153930, 6.13838153640245e-7, -33.8074190956349, -0.0789797357153930, -6.13838153640245e-7, -29.2315964718309]
  1: [4.76837158203125e-7, -2.88366995014816, 0.000000000000000, -4.76837158203125e-7, 2.88366995014817, 14.4183480425128]
  2: [0.000000000000000, -2.96264968586362, 14.8132484293181, 0.000000000000000, 2.96264968586362, 0.000000000000000]

print_indexed (map2 (chop(), map2 (lambda x, y : abs (chop()(x) - chop()(y))/abs (chop()(x)) * 100, ef, ef2)))

  0: [0.00616427213572451, 6.13838154611604e-7, 14.0940184048878,
0.00616427213572451, 6.13838154611604e-7, 18.2549943272073]
  1: [4.76837155929388e-7, 5.16212444975993, 0.000000000000000,
4.76837155929388e-7, 2.00063237766665, 6.53331518958011]
  2: [0.000000000000000, 2.16075060864196, 24.4603272876847, 0.000000000000000,
2.62441648235529, 0.000000000000000]

  0: [0.00616427213572451, 6.13838154611604e-7, 14.0940184048878, 0.00616427213572451, 6.13838154611604e-7, 18.2549943272073]
  1: [4.76837155929388e-7, 5.16212444975993, 0.000000000000000, 4.76837155929388e-7, 2.00063237766665, 6.53331518958011]
  2: [0.000000000000000, 2.16075060864196, 24.4603272876847, 0.000000000000000, 2.62441648235529, 0.000000000000000]

nds_xy = make_nodes (17, (0., 0.), (0., 3.25))
nds_xy += make_nodes (17, (5.5, 0.), (0., 3.25))
print 'nodes:'
print_indexed (nds_xy, '  ')
print 

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 'elements:'
print_indexed (els, '  ')
print

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))]
print 'supports:', supps
print

ch_stup1d =  (3e7, 1000, 0.6*0.8^3/12)
ch_stup1g =  (3e7, 1000, 0.4*0.6^3/12)
ch_stup2d =  (3e7, 1000, 0.6*0.8^3/12)
ch_stup2g =  (3e7, 1000, 0.4*0.6^3/12)
ch_greda =  (3e7, 1000, 0.4*0.8^3/12)

chars = [ch_stup1d, ch_stup1g, ch_stup2d, ch_stup2g, ch_greda]
print 'characteristics:'
print_list (chars, '   ')
print

echs = [0 for i in stup1d_e] + [1 for i in stup1g_e]
echs += [2 for i in stup2d_e] + [3 for i in stup2g_e]
echs += [4 for i in grede_e]
print 'elements <-> chars:', echs
print 

ndlds = [(i, (100./16*i, 0., 0.)) for i in stup1_n[1:]]
print 'nodal loads:'
print_list (ndlds, '   ')
print

ellds = [(i, df2, (0., 50.)) for i in grede_e]
print 'element loads:'
print_list (ellds, '   ')

WARNING: Output truncated!

full_output.txt




nodes:
    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]

elements:
    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]

...

   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]

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

characteristics:
   (3.00000000000000e7, 1000, 0.0256000000000000)
   (3.00000000000000e7, 1000, 0.00720000000000000)
   (3.00000000000000e7, 1000, 0.0256000000000000)
   (3.00000000000000e7, 1000, 0.00720000000000000)
   (3.00000000000000e7, 1000, 0.0170666666666667)

elements <-> chars: [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,
2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4]

nodal loads:
   (1, (6.25000000000000, 0.000000000000000, 0.000000000000000))
   (2, (12.5000000000000, 0.000000000000000, 0.000000000000000))
   (3, (18.7500000000000, 0.000000000000000, 0.000000000000000))
   (4, (25.0000000000000, 0.000000000000000, 0.000000000000000))
   (5, (31.2500000000000, 0.000000000000000, 0.000000000000000))
   (6, (37.5000000000000, 0.000000000000000, 0.000000000000000))
   (7, (43.7500000000000, 0.000000000000000, 0.000000000000000))
   (8, (50.0000000000000, 0.000000000000000, 0.000000000000000))
   (9, (56.2500000000000, 0.000000000000000, 0.000000000000000))
   (10, (62.5000000000000, 0.000000000000000, 0.000000000000000))
   (11, (68.7500000000000, 0.000000000000000, 0.000000000000000))
   (12, (75.0000000000000, 0.000000000000000, 0.000000000000000))
   (13, (81.2500000000000, 0.000000000000000, 0.000000000000000))
   (14, (87.5000000000000, 0.000000000000000, 0.000000000000000))
   (15, (93.7500000000000, 0.000000000000000, 0.000000000000000))
   (16, (100.000000000000, 0.000000000000000, 0.000000000000000))

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

WARNING: Output truncated!

full_output.txt




nodes:
    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]

elements:
    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]

...

   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]

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

characteristics:
   (3.00000000000000e7, 1000, 0.0256000000000000)
   (3.00000000000000e7, 1000, 0.00720000000000000)
   (3.00000000000000e7, 1000, 0.0256000000000000)
   (3.00000000000000e7, 1000, 0.00720000000000000)
   (3.00000000000000e7, 1000, 0.0170666666666667)

elements <-> chars: [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]

nodal loads:
   (1, (6.25000000000000, 0.000000000000000, 0.000000000000000))
   (2, (12.5000000000000, 0.000000000000000, 0.000000000000000))
   (3, (18.7500000000000, 0.000000000000000, 0.000000000000000))
   (4, (25.0000000000000, 0.000000000000000, 0.000000000000000))
   (5, (31.2500000000000, 0.000000000000000, 0.000000000000000))
   (6, (37.5000000000000, 0.000000000000000, 0.000000000000000))
   (7, (43.7500000000000, 0.000000000000000, 0.000000000000000))
   (8, (50.0000000000000, 0.000000000000000, 0.000000000000000))
   (9, (56.2500000000000, 0.000000000000000, 0.000000000000000))
   (10, (62.5000000000000, 0.000000000000000, 0.000000000000000))
   (11, (68.7500000000000, 0.000000000000000, 0.000000000000000))
   (12, (75.0000000000000, 0.000000000000000, 0.000000000000000))
   (13, (81.2500000000000, 0.000000000000000, 0.000000000000000))
   (14, (87.5000000000000, 0.000000000000000, 0.000000000000000))
   (15, (93.7500000000000, 0.000000000000000, 0.000000000000000))
   (16, (100.000000000000, 0.000000000000000, 0.000000000000000))

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

full_output.txt

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds, cf_scale = 0.25, df_scale = 0.02)

nds = coords_left_to_right (nds_xy)

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

K = structure_stiffness_matrix (f2k, f2r, ndof, tbl_dof, nds, els, echs, chars)
q = structure_force_vector (ndof, tbl_dof, ndlds, ellds, els, nds, f2r)

u = nodal_displacements (K, q, tbl_dof)
print 'nodal displacements:'
print_indexed (map2 (ND(6), u), '  ')

nodal displacements:
    0: [0.000000, 0.000000, 0.000000]
    1: [0.00460252, -3.15704e-7, -0.00190989]
    2: [0.0129051, -6.05424e-7, -0.00232340]
    3: [0.0219816, -8.59790e-7, -0.00239105]
    4: [0.0310687, -1.07738e-6, -0.00235369]
    5: [0.0399279, -1.25902e-6, -0.00227861]
    6: [0.0484498, -1.40634e-6, -0.00218338]
    7: [0.0565632, -1.52148e-6, -0.00206495]
    8: [0.0642335, -1.60690e-6, -0.00198240]
    9: [0.0743132, -1.66505e-6, -0.00187855]
   10: [0.0834441, -1.69987e-6, -0.00168890]
   11: [0.0915660, -1.71522e-6, -0.00150053]
   12: [0.0985823, -1.71530e-6, -0.00129109]
   13: [0.104393, -1.70473e-6, -0.00106293]
   14: [0.108898, -1.68849e-6, -0.000822982]
   15: [0.112001, -1.67198e-6, -0.000534742]
   16: [0.113648, -1.66095e-6, -0.000424804]
   17: [0.000000, 0.000000, 0.000000]
   18: [0.00460253, 7.92371e-7, -0.00180638]
   19: [0.0129051, 1.52897e-6, -0.00224481]
   20: [0.0219816, 2.20041e-6, -0.00230650]
   21: [0.0310687, 2.80530e-6, -0.00227040]
   22: [0.0399279, 3.34443e-6, -0.00219576]
   23: [0.0484498, 3.81946e-6, -0.00209732]
   24: [0.0565632, 4.23252e-6, -0.00199256]
   25: [0.0642335, 4.58607e-6, -0.00185317]
   26: [0.0743132, 4.88255e-6, -0.00159791]
   27: [0.0834441, 5.12591e-6, -0.00143732]
   28: [0.0915660, 5.32001e-6, -0.00124344]
   29: [0.0985823, 5.46905e-6, -0.00103468]
   30: [0.104393, 5.57764e-6, -0.000808337]
   31: [0.108898, 5.65078e-6, -0.000557903]
   32: [0.112001, 5.69386e-6, -0.000324483]
   33: [0.113648, 5.71261e-6, 0.0000710427]

nodal displacements:
    0: [0.000000, 0.000000, 0.000000]
    1: [0.00460252, -3.15704e-7, -0.00190989]
    2: [0.0129051, -6.05424e-7, -0.00232340]
    3: [0.0219816, -8.59790e-7, -0.00239105]
    4: [0.0310687, -1.07738e-6, -0.00235369]
    5: [0.0399279, -1.25902e-6, -0.00227861]
    6: [0.0484498, -1.40634e-6, -0.00218338]
    7: [0.0565632, -1.52148e-6, -0.00206495]
    8: [0.0642335, -1.60690e-6, -0.00198240]
    9: [0.0743132, -1.66505e-6, -0.00187855]
   10: [0.0834441, -1.69987e-6, -0.00168890]
   11: [0.0915660, -1.71522e-6, -0.00150053]
   12: [0.0985823, -1.71530e-6, -0.00129109]
   13: [0.104393, -1.70473e-6, -0.00106293]
   14: [0.108898, -1.68849e-6, -0.000822982]
   15: [0.112001, -1.67198e-6, -0.000534742]
   16: [0.113648, -1.66095e-6, -0.000424804]
   17: [0.000000, 0.000000, 0.000000]
   18: [0.00460253, 7.92371e-7, -0.00180638]
   19: [0.0129051, 1.52897e-6, -0.00224481]
   20: [0.0219816, 2.20041e-6, -0.00230650]
   21: [0.0310687, 2.80530e-6, -0.00227040]
   22: [0.0399279, 3.34443e-6, -0.00219576]
   23: [0.0484498, 3.81946e-6, -0.00209732]
   24: [0.0565632, 4.23252e-6, -0.00199256]
   25: [0.0642335, 4.58607e-6, -0.00185317]
   26: [0.0743132, 4.88255e-6, -0.00159791]
   27: [0.0834441, 5.12591e-6, -0.00143732]
   28: [0.0915660, 5.32001e-6, -0.00124344]
   29: [0.0985823, 5.46905e-6, -0.00103468]
   30: [0.104393, 5.57764e-6, -0.000808337]
   31: [0.108898, 5.65078e-6, -0.000557903]
   32: [0.112001, 5.69386e-6, -0.000324483]
   33: [0.113648, 5.71261e-6, 0.0000710427]

uu = extract_from (u, ix = 0)
i_mu, u_max = iamax (uu)
i_mu, u_max

(16, 0.113648441675127)

(16, 0.113648441675127)

ef = all_end_forces (f2k, f2r, u, nds, els, echs, chars, ellds)
print 'element end forces:'
print_indexed (map2 (ND(6), ef), '  ')

element end forces:
    0: [-2914.19, -402.422, 1105.26, 2914.19, 402.422, 202.616]
    1: [-2674.33, -382.155, 718.718, 2674.33, 382.155, 523.285]
    2: [-2347.99, -380.039, 633.549, 2347.99, 380.039, 601.578]
    3: [-2008.55, -369.639, 591.837, 2008.55, 369.639, 609.491]
    4: [-1676.60, -357.510, 563.211, 1676.60, 357.510, 598.696]
    5: [-1359.88, -341.281, 532.078, 1359.88, 341.281, 577.087]
    6: [-1062.86, -324.811, 499.834, 1062.86, 324.811, 555.803]
    7: [-788.513, -293.521, 457.463, 788.513, 293.521, 496.479]
    8: [-536.726, -287.355, 460.050, 536.726, 287.355, 473.853]
    9: [-321.420, -251.724, 396.446, 321.420, 251.724, 421.655]
   10: [-141.693, -221.918, 348.098, 141.693, 221.918, 373.136]
   11: [-0.792179, -187.247, 290.357, 0.792179, 187.247, 318.197]
   12: [97.6329, -149.901, 228.425, -97.6329, 149.901, 258.752]
   13: [149.898, -108.744, 160.761, -149.898, 108.744, 192.656]
   14: [152.361, -67.7133, 90.8772, -152.361, 67.7133, 129.191]
   15: [101.847, -6.68146, 3.55068, -101.847, 6.68146, 18.1641]
   16: [7314.19, -447.578, 1154.18, -7314.19, 447.578, 300.451]
   17: [6799.33, -461.595, 853.695, -6799.33, 461.595, 646.489]
   18: [6197.99, -451.211, 747.796, -6197.99, 451.211, 718.639]
   19: [5583.55, -442.861, 711.118, -5583.55, 442.861, 728.180]
   20: [4976.60, -429.990, 681.097, -4976.60, 429.990, 716.372]
   21: [4384.88, -414.969, 651.061, -4384.88, 414.969, 697.587]
   22: [3812.86, -393.939, 615.395, -3812.86, 393.939, 664.906]
   23: [3263.51, -381.479, 586.965, -3263.51, 381.479, 652.843]
   24: [2736.73, -337.645, 531.709, -2736.73, 337.645, 565.638]
   25: [2246.42, -317.026, 504.495, -2246.42, 317.026, 525.841]
   26: [1791.69, -284.332, 449.154, -1791.69, 284.332, 474.925]
   27: [1375.79, -250.253, 392.786, -1375.79, 250.253, 420.535]
   28: [1002.37, -212.599, 330.431, -1002.37, 212.599, 360.517]
   29: [675.102, -172.506, 263.679, -675.102, 172.506, 296.967]
   30: [397.639, -126.037, 189.296, -397.639, 126.037, 220.323]
   31: [173.153, -93.3185, 125.355, -173.153, 93.3185, 177.930]
   32: [-14.0176, 239.860, -921.334, 14.0176, -514.860, -1154.15]
   33: [10.3844, 326.340, -1156.83, -10.3844, -601.340, -1394.28]
   34: [8.34988, 339.440, -1193.42, -8.34988, -614.440, -1429.76]
   35: [12.8705, 331.951, -1172.70, -12.8705, -606.951, -1409.28]
   36: [15.0219, 316.719, -1130.77, -15.0219, -591.719, -1367.43]
   37: [21.0300, 297.028, -1076.92, -21.0300, -572.028, -1312.98]
   38: [12.4591, 274.343, -1013.27, -12.4591, -549.343, -1251.87]
   39: [43.8342, 251.787, -956.529, -43.8342, -526.787, -1184.55]
   40: [20.6188, 215.306, -870.300, -20.6188, -490.306, -1070.13]
   41: [32.6946, 179.727, -769.753, -32.6946, -454.727, -974.995]
   42: [34.0792, 140.901, -663.493, -34.0792, -415.901, -867.711]
   43: [37.6533, 98.4250, -546.621, -37.6533, -373.425, -750.966]
   44: [40.0931, 52.2653, -419.513, -40.0931, -327.265, -624.196]
   45: [46.4697, 2.46296, -283.533, -46.4697, -277.463, -486.263]
   46: [32.7182, -50.5145, -132.742, -32.7182, -224.485, -345.678]
   47: [93.3185, -101.847, -18.1641, -93.3185, -173.153, -177.930]

element end forces:
    0: [-2914.19, -402.422, 1105.26, 2914.19, 402.422, 202.616]
    1: [-2674.33, -382.155, 718.718, 2674.33, 382.155, 523.285]
    2: [-2347.99, -380.039, 633.549, 2347.99, 380.039, 601.578]
    3: [-2008.55, -369.639, 591.837, 2008.55, 369.639, 609.491]
    4: [-1676.60, -357.510, 563.211, 1676.60, 357.510, 598.696]
    5: [-1359.88, -341.281, 532.078, 1359.88, 341.281, 577.087]
    6: [-1062.86, -324.811, 499.834, 1062.86, 324.811, 555.803]
    7: [-788.513, -293.521, 457.463, 788.513, 293.521, 496.479]
    8: [-536.726, -287.355, 460.050, 536.726, 287.355, 473.853]
    9: [-321.420, -251.724, 396.446, 321.420, 251.724, 421.655]
   10: [-141.693, -221.918, 348.098, 141.693, 221.918, 373.136]
   11: [-0.792179, -187.247, 290.357, 0.792179, 187.247, 318.197]
   12: [97.6329, -149.901, 228.425, -97.6329, 149.901, 258.752]
   13: [149.898, -108.744, 160.761, -149.898, 108.744, 192.656]
   14: [152.361, -67.7133, 90.8772, -152.361, 67.7133, 129.191]
   15: [101.847, -6.68146, 3.55068, -101.847, 6.68146, 18.1641]
   16: [7314.19, -447.578, 1154.18, -7314.19, 447.578, 300.451]
   17: [6799.33, -461.595, 853.695, -6799.33, 461.595, 646.489]
   18: [6197.99, -451.211, 747.796, -6197.99, 451.211, 718.639]
   19: [5583.55, -442.861, 711.118, -5583.55, 442.861, 728.180]
   20: [4976.60, -429.990, 681.097, -4976.60, 429.990, 716.372]
   21: [4384.88, -414.969, 651.061, -4384.88, 414.969, 697.587]
   22: [3812.86, -393.939, 615.395, -3812.86, 393.939, 664.906]
   23: [3263.51, -381.479, 586.965, -3263.51, 381.479, 652.843]
   24: [2736.73, -337.645, 531.709, -2736.73, 337.645, 565.638]
   25: [2246.42, -317.026, 504.495, -2246.42, 317.026, 525.841]
   26: [1791.69, -284.332, 449.154, -1791.69, 284.332, 474.925]
   27: [1375.79, -250.253, 392.786, -1375.79, 250.253, 420.535]
   28: [1002.37, -212.599, 330.431, -1002.37, 212.599, 360.517]
   29: [675.102, -172.506, 263.679, -675.102, 172.506, 296.967]
   30: [397.639, -126.037, 189.296, -397.639, 126.037, 220.323]
   31: [173.153, -93.3185, 125.355, -173.153, 93.3185, 177.930]
   32: [-14.0176, 239.860, -921.334, 14.0176, -514.860, -1154.15]
   33: [10.3844, 326.340, -1156.83, -10.3844, -601.340, -1394.28]
   34: [8.34988, 339.440, -1193.42, -8.34988, -614.440, -1429.76]
   35: [12.8705, 331.951, -1172.70, -12.8705, -606.951, -1409.28]
   36: [15.0219, 316.719, -1130.77, -15.0219, -591.719, -1367.43]
   37: [21.0300, 297.028, -1076.92, -21.0300, -572.028, -1312.98]
   38: [12.4591, 274.343, -1013.27, -12.4591, -549.343, -1251.87]
   39: [43.8342, 251.787, -956.529, -43.8342, -526.787, -1184.55]
   40: [20.6188, 215.306, -870.300, -20.6188, -490.306, -1070.13]
   41: [32.6946, 179.727, -769.753, -32.6946, -454.727, -974.995]
   42: [34.0792, 140.901, -663.493, -34.0792, -415.901, -867.711]
   43: [37.6533, 98.4250, -546.621, -37.6533, -373.425, -750.966]
   44: [40.0931, 52.2653, -419.513, -40.0931, -327.265, -624.196]
   45: [46.4697, 2.46296, -283.533, -46.4697, -277.463, -486.263]
   46: [32.7182, -50.5145, -132.742, -32.7182, -224.485, -345.678]
   47: [93.3185, -101.847, -18.1641, -93.3185, -173.153, -177.930]

Hs = extract_from (ef, ix = 3)
eHs = range (len (Hs))

i_mn, N_max = iamax (Hs)
i_mn, N_max

(16, -7314.19398515837)

(16, -7314.19398515837)

Kg = structure_stiffness_matrix (f2kg, f2r, ndof, tbl_dof, nds, els, eHs, Hs)
K2 = K + Kg

elldsg = [(i, df2g, (50., Hs[i], ch_greda[0], ch_greda[2])) for i in grede_e]
qg = structure_force_vector (ndof, tbl_dof, [], elldsg, els, nds, f2r)
q2 = q + qg

u2 = nodal_displacements (K2, q2, tbl_dof)
print 'nodal displacements:'
print_indexed (map2 (ND(6), u2), '  ')

nodal displacements:
    0: [0.000000, 0.000000, 0.000000]
    1: [0.00464879, -3.21191e-7, -0.00193030]
    2: [0.0130525, -6.15941e-7, -0.00235259]
    3: [0.0222448, -8.74706e-7, -0.00242097]
    4: [0.0314444, -1.09605e-6, -0.00238172]
    5: [0.0404065, -1.28082e-6, -0.00230381]
    6: [0.0490202, -1.43074e-6, -0.00220564]
    7: [0.0572144, -1.54799e-6, -0.00208459]
    8: [0.0649585, -1.63510e-6, -0.00200179]
    9: [0.0751475, -1.69452e-6, -0.00189777]
   10: [0.0843696, -1.73022e-6, -0.00170383]
   11: [0.0925615, -1.74615e-6, -0.00151165]
   12: [0.0996284, -1.74658e-6, -0.00129880]
   13: [0.105473, -1.73618e-6, -0.00106780]
   14: [0.109998, -1.72003e-6, -0.000825641]
   15: [0.113110, -1.70355e-6, -0.000535825]
   16: [0.114761, -1.69253e-6, -0.000425213]
   17: [0.000000, 0.000000, 0.000000]
   18: [0.00464879, 7.97858e-7, -0.00182738]
   19: [0.0130525, 1.53948e-6, -0.00227308]
   20: [0.0222448, 2.21533e-6, -0.00233565]
   21: [0.0314444, 2.82397e-6, -0.00229768]
   22: [0.0404065, 3.36624e-6, -0.00222035]
   23: [0.0490202, 3.84386e-6, -0.00211902]
   24: [0.0572144, 4.25904e-6, -0.00201195]
   25: [0.0649585, 4.61427e-6, -0.00187155]
   26: [0.0751475, 4.91202e-6, -0.00161401]
   27: [0.0843696, 5.15626e-6, -0.00145052]
   28: [0.0925615, 5.35094e-6, -0.00125331]
   29: [0.0996284, 5.50033e-6, -0.00104163]
   30: [0.105473, 5.60910e-6, -0.000812778]
   31: [0.109998, 5.68232e-6, -0.000560323]
   32: [0.113110, 5.72543e-6, -0.000325530]
   33: [0.114761, 5.74419e-6, 0.0000708476]

nodal displacements:
    0: [0.000000, 0.000000, 0.000000]
    1: [0.00464879, -3.21191e-7, -0.00193030]
    2: [0.0130525, -6.15941e-7, -0.00235259]
    3: [0.0222448, -8.74706e-7, -0.00242097]
    4: [0.0314444, -1.09605e-6, -0.00238172]
    5: [0.0404065, -1.28082e-6, -0.00230381]
    6: [0.0490202, -1.43074e-6, -0.00220564]
    7: [0.0572144, -1.54799e-6, -0.00208459]
    8: [0.0649585, -1.63510e-6, -0.00200179]
    9: [0.0751475, -1.69452e-6, -0.00189777]
   10: [0.0843696, -1.73022e-6, -0.00170383]
   11: [0.0925615, -1.74615e-6, -0.00151165]
   12: [0.0996284, -1.74658e-6, -0.00129880]
   13: [0.105473, -1.73618e-6, -0.00106780]
   14: [0.109998, -1.72003e-6, -0.000825641]
   15: [0.113110, -1.70355e-6, -0.000535825]
   16: [0.114761, -1.69253e-6, -0.000425213]
   17: [0.000000, 0.000000, 0.000000]
   18: [0.00464879, 7.97858e-7, -0.00182738]
   19: [0.0130525, 1.53948e-6, -0.00227308]
   20: [0.0222448, 2.21533e-6, -0.00233565]
   21: [0.0314444, 2.82397e-6, -0.00229768]
   22: [0.0404065, 3.36624e-6, -0.00222035]
   23: [0.0490202, 3.84386e-6, -0.00211902]
   24: [0.0572144, 4.25904e-6, -0.00201195]
   25: [0.0649585, 4.61427e-6, -0.00187155]
   26: [0.0751475, 4.91202e-6, -0.00161401]
   27: [0.0843696, 5.15626e-6, -0.00145052]
   28: [0.0925615, 5.35094e-6, -0.00125331]
   29: [0.0996284, 5.50033e-6, -0.00104163]
   30: [0.105473, 5.60910e-6, -0.000812778]
   31: [0.109998, 5.68232e-6, -0.000560323]
   32: [0.113110, 5.72543e-6, -0.000325530]
   33: [0.114761, 5.74419e-6, 0.0000708476]

uu2 = extract_from (u2, ix = 0)
i_mu2, u2_max = iamax (uu2)
i_mu2, u2_max

(16, 0.114761314293619)

(16, 0.114761314293619)

abs (u_max - u2_max) / abs (u_max) * 100

0.979223825764268

0.979223825764268

ef2e = all_end_forces (f2k, f2r, u2, nds, els, echs, chars, ellds)
ef2g = all_end_forces (f2kg, f2r, u2, nds, els, eHs, Hs, elldsg)
ef2 = map (lambda x, y : x + y, ef2e, ef2g)
print 'element end forces:'
print_indexed (map2 (ND(6), map2 (chop(), ef2)), '  ')

element end forces:
    0: [-2964.84, -410.379, 1117.76, 2964.84, 410.379, 202.425]
    1: [-2720.77, -394.814, 730.433, 2720.77, 394.814, 530.238]
    2: [-2388.60, -392.186, 642.712, 2388.60, 392.186, 610.309]
    3: [-2043.17, -380.428, 599.660, 2043.17, 380.428, 618.254]
    4: [-1705.61, -366.706, 569.939, 1705.61, 366.706, 606.831]
    5: [-1383.84, -348.903, 537.877, 1383.84, 348.903, 584.345]
    6: [-1082.36, -330.995, 504.872, 1082.36, 330.995, 562.154]
    7: [-804.059, -298.243, 462.010, 804.059, 298.243, 501.175]
    8: [-548.455, -292.677, 465.937, 548.455, 292.677, 479.793]
    9: [-329.572, -255.396, 400.630, 329.572, 255.396, 426.444]
   10: [-147.000, -224.395, 351.283, 147.000, 224.395, 376.842]
   11: [-3.97948, -188.762, 292.590, 3.97948, 188.762, 320.883]
   12: [95.9233, -150.723, 229.863, -95.9233, 150.723, 260.557]
   13: [149.124, -109.110, 161.558, -149.124, 109.110, 193.726]
   14: [152.094, -67.8302, 91.2117, -152.094, 67.8302, 129.711]
   15: [101.791, -6.69056, 3.60781, -101.791, 6.69056, 18.3046]
   16: [7364.84, -439.621, 1159.59, -7364.84, 439.621, 303.182]
   17: [6845.77, -448.936, 862.593, -6845.77, 448.936, 653.588]
   18: [6238.60, -439.064, 756.646, -6238.60, 439.064, 727.286]
   19: [5618.17, -432.072, 718.885, -5618.17, 432.072, 736.715]
   20: [5005.61, -420.794, 687.920, -5005.61, 420.794, 724.260]
   21: [4408.84, -407.347, 656.998, -4408.84, 407.347, 704.648]
   22: [3832.36, -387.755, 620.531, -3832.36, 387.755, 670.915]
   23: [3279.06, -376.757, 591.812, -3279.06, 376.757, 657.919]
   24: [2748.46, -332.323, 537.042, -2748.46, 332.323, 570.893]
   25: [2254.57, -313.354, 508.792, -2254.57, 313.354, 530.324]
   26: [1797.00, -281.855, 452.341, -1797.00, 281.855, 478.364]
   27: [1378.98, -248.738, 395.070, -1378.98, 248.738, 423.050]
   28: [1004.08, -211.777, 331.919, -1004.08, 211.777, 362.215]
   29: [675.876, -172.140, 264.523, -675.876, 172.140, 297.988]
   30: [397.906, -125.920, 189.659, -397.906, 125.920, 220.818]
   31: [173.209, -93.3094, 125.445, -173.209, 93.3094, 178.096]
   32: [-9.31502, 244.070, -932.858, 9.31502, -519.070, -1165.78]
   33: [9.87202, 332.170, -1172.95, -9.87202, -607.170, -1410.23]
   34: [6.99214, 345.434, -1209.97, -6.99214, -620.434, -1446.17]
   35: [11.2782, 337.560, -1188.19, -11.2782, -612.560, -1424.63]
   36: [13.4469, 321.767, -1144.71, -13.4469, -596.767, -1381.26]
   37: [19.5921, 301.481, -1089.22, -19.5921, -576.481, -1325.18]
   38: [10.9981, 278.298, -1024.16, -10.9981, -553.298, -1262.73]
   39: [44.4331, 255.604, -967.112, -44.4331, -530.604, -1194.96]
   40: [18.9698, 218.883, -880.423, -18.9698, -493.883, -1079.69]
   41: [31.4991, 182.571, -777.727, -31.4991, -457.571, -982.665]
   42: [33.1171, 143.021, -669.431, -33.1171, -418.021, -873.434]
   43: [36.9605, 99.9028, -550.746, -36.9605, -374.903, -754.969]
   44: [39.6367, 53.2005, -422.115, -39.6367, -328.201, -626.738]
   45: [46.2207, 2.97006, -284.938, -46.2207, -277.970, -487.647]
   46: [32.6103, -50.3032, -133.319, -32.6103, -224.697, -346.263]
   47: [93.3094, -101.791, -18.3046, -93.3094, -173.209, -178.096]

element end forces:
    0: [-2964.84, -410.379, 1117.76, 2964.84, 410.379, 202.425]
    1: [-2720.77, -394.814, 730.433, 2720.77, 394.814, 530.238]
    2: [-2388.60, -392.186, 642.712, 2388.60, 392.186, 610.309]
    3: [-2043.17, -380.428, 599.660, 2043.17, 380.428, 618.254]
    4: [-1705.61, -366.706, 569.939, 1705.61, 366.706, 606.831]
    5: [-1383.84, -348.903, 537.877, 1383.84, 348.903, 584.345]
    6: [-1082.36, -330.995, 504.872, 1082.36, 330.995, 562.154]
    7: [-804.059, -298.243, 462.010, 804.059, 298.243, 501.175]
    8: [-548.455, -292.677, 465.937, 548.455, 292.677, 479.793]
    9: [-329.572, -255.396, 400.630, 329.572, 255.396, 426.444]
   10: [-147.000, -224.395, 351.283, 147.000, 224.395, 376.842]
   11: [-3.97948, -188.762, 292.590, 3.97948, 188.762, 320.883]
   12: [95.9233, -150.723, 229.863, -95.9233, 150.723, 260.557]
   13: [149.124, -109.110, 161.558, -149.124, 109.110, 193.726]
   14: [152.094, -67.8302, 91.2117, -152.094, 67.8302, 129.711]
   15: [101.791, -6.69056, 3.60781, -101.791, 6.69056, 18.3046]
   16: [7364.84, -439.621, 1159.59, -7364.84, 439.621, 303.182]
   17: [6845.77, -448.936, 862.593, -6845.77, 448.936, 653.588]
   18: [6238.60, -439.064, 756.646, -6238.60, 439.064, 727.286]
   19: [5618.17, -432.072, 718.885, -5618.17, 432.072, 736.715]
   20: [5005.61, -420.794, 687.920, -5005.61, 420.794, 724.260]
   21: [4408.84, -407.347, 656.998, -4408.84, 407.347, 704.648]
   22: [3832.36, -387.755, 620.531, -3832.36, 387.755, 670.915]
   23: [3279.06, -376.757, 591.812, -3279.06, 376.757, 657.919]
   24: [2748.46, -332.323, 537.042, -2748.46, 332.323, 570.893]
   25: [2254.57, -313.354, 508.792, -2254.57, 313.354, 530.324]
   26: [1797.00, -281.855, 452.341, -1797.00, 281.855, 478.364]
   27: [1378.98, -248.738, 395.070, -1378.98, 248.738, 423.050]
   28: [1004.08, -211.777, 331.919, -1004.08, 211.777, 362.215]
   29: [675.876, -172.140, 264.523, -675.876, 172.140, 297.988]
   30: [397.906, -125.920, 189.659, -397.906, 125.920, 220.818]
   31: [173.209, -93.3094, 125.445, -173.209, 93.3094, 178.096]
   32: [-9.31502, 244.070, -932.858, 9.31502, -519.070, -1165.78]
   33: [9.87202, 332.170, -1172.95, -9.87202, -607.170, -1410.23]
   34: [6.99214, 345.434, -1209.97, -6.99214, -620.434, -1446.17]
   35: [11.2782, 337.560, -1188.19, -11.2782, -612.560, -1424.63]
   36: [13.4469, 321.767, -1144.71, -13.4469, -596.767, -1381.26]
   37: [19.5921, 301.481, -1089.22, -19.5921, -576.481, -1325.18]
   38: [10.9981, 278.298, -1024.16, -10.9981, -553.298, -1262.73]
   39: [44.4331, 255.604, -967.112, -44.4331, -530.604, -1194.96]
   40: [18.9698, 218.883, -880.423, -18.9698, -493.883, -1079.69]
   41: [31.4991, 182.571, -777.727, -31.4991, -457.571, -982.665]
   42: [33.1171, 143.021, -669.431, -33.1171, -418.021, -873.434]
   43: [36.9605, 99.9028, -550.746, -36.9605, -374.903, -754.969]
   44: [39.6367, 53.2005, -422.115, -39.6367, -328.201, -626.738]
   45: [46.2207, 2.97006, -284.938, -46.2207, -277.970, -487.647]
   46: [32.6103, -50.3032, -133.319, -32.6103, -224.697, -346.263]
   47: [93.3094, -101.791, -18.3046, -93.3094, -173.209, -178.096]

print_indexed (zip (extract_from (ef, ix=3), extract_from (ef2, ix=3)))

  0: (2914.19398515826, 2964.83921638775)
  1: (2674.33389991300, 2720.76956173476)
  2: (2347.99410437279, 2388.59959808417)
  3: (2008.55368531835, 2043.16519301094)
  4: (1676.60316749151, 1705.60558942421)
  5: (1359.88385255614, 1383.83897325495)
  6: (1062.85616430883, 1082.35772757289)
  7: (788.513165393053, 804.059342316739)
  8: (536.725754133138, 548.455111013631)
  9: (321.419835254708, 329.571693870204)
 10: (141.692892440267, 147.000410656805)
 11: (0.792178794077699, 3.97948000324686)
 12: (-97.6328601208534, -95.9232972594873)
 13: (-149.898118690480, -149.123803489092)
 14: (-152.361078772023, -152.093867236119)
 15: (-101.846560030863, -101.790650151219)
 16: (-7314.19398515837, -7364.83921638747)
 17: (-6799.33389991310, -6845.76956173448)
 18: (-6197.99410437291, -6238.59959808390)
 19: (-5583.55368531845, -5618.16519301066)
 20: (-4976.60316749160, -5005.60558942390)
 21: (-4384.88385255622, -4408.83897325466)
 22: (-3812.85616430888, -3832.35772757261)
 23: (-3263.51316539308, -3279.05934231640)
 24: (-2736.72575413316, -2748.45511101332)
 25: (-2246.41983525470, -2254.57169386991)
 26: (-1791.69289244027, -1797.00041065663)
 27: (-1375.79217879406, -1378.97948000312)
 28: (-1002.36713987907, -1004.07670274046)
 29: (-675.101881309492, -675.876196510864)
 30: (-397.638921227961, -397.906132763870)
 31: (-173.153439969137, -173.209349848752)
 32: (14.0175611339509, 9.31502197310328)
 33: (-10.3843817263842, -9.87202137708664)
 34: (-8.34987826645374, -6.99213789403439)
 35: (-12.8705210387707, -11.2782401442528)
 36: (-15.0218656361103, -13.4468937814236)
 37: (-21.0299560725689, -19.5920789539814)
 38: (-12.4591064453125, -10.9980939626694)
 39: (-43.8341569900513, -44.4330536723137)
 40: (-20.6188215017319, -18.9697536230087)
 41: (-32.6945824027061, -31.4991106390953)
 42: (-34.0791532993317, -33.1170613765717)
 43: (-37.6532667875290, -36.9604554176331)
 44: (-40.0930626392365, -39.6366543769836)
 45: (-46.4696916341782, -46.2206619977951)
 46: (-32.7181651592255, -32.6103112697601)
 47: (-93.3185379505157, -93.3094445466995)

  0: (2914.19398515826, 2964.83921638775)
  1: (2674.33389991300, 2720.76956173476)
  2: (2347.99410437279, 2388.59959808417)
  3: (2008.55368531835, 2043.16519301094)
  4: (1676.60316749151, 1705.60558942421)
  5: (1359.88385255614, 1383.83897325495)
  6: (1062.85616430883, 1082.35772757289)
  7: (788.513165393053, 804.059342316739)
  8: (536.725754133138, 548.455111013631)
  9: (321.419835254708, 329.571693870204)
 10: (141.692892440267, 147.000410656805)
 11: (0.792178794077699, 3.97948000324686)
 12: (-97.6328601208534, -95.9232972594873)
 13: (-149.898118690480, -149.123803489092)
 14: (-152.361078772023, -152.093867236119)
 15: (-101.846560030863, -101.790650151219)
 16: (-7314.19398515837, -7364.83921638747)
 17: (-6799.33389991310, -6845.76956173448)
 18: (-6197.99410437291, -6238.59959808390)
 19: (-5583.55368531845, -5618.16519301066)
 20: (-4976.60316749160, -5005.60558942390)
 21: (-4384.88385255622, -4408.83897325466)
 22: (-3812.85616430888, -3832.35772757261)
 23: (-3263.51316539308, -3279.05934231640)
 24: (-2736.72575413316, -2748.45511101332)
 25: (-2246.41983525470, -2254.57169386991)
 26: (-1791.69289244027, -1797.00041065663)
 27: (-1375.79217879406, -1378.97948000312)
 28: (-1002.36713987907, -1004.07670274046)
 29: (-675.101881309492, -675.876196510864)
 30: (-397.638921227961, -397.906132763870)
 31: (-173.153439969137, -173.209349848752)
 32: (14.0175611339509, 9.31502197310328)
 33: (-10.3843817263842, -9.87202137708664)
 34: (-8.34987826645374, -6.99213789403439)
 35: (-12.8705210387707, -11.2782401442528)
 36: (-15.0218656361103, -13.4468937814236)
 37: (-21.0299560725689, -19.5920789539814)
 38: (-12.4591064453125, -10.9980939626694)
 39: (-43.8341569900513, -44.4330536723137)
 40: (-20.6188215017319, -18.9697536230087)
 41: (-32.6945824027061, -31.4991106390953)
 42: (-34.0791532993317, -33.1170613765717)
 43: (-37.6532667875290, -36.9604554176331)
 44: (-40.0930626392365, -39.6366543769836)
 45: (-46.4696916341782, -46.2206619977951)
 46: (-32.7181651592255, -32.6103112697601)
 47: (-93.3185379505157, -93.3094445466995)

Hs2 = extract_from (ef2, ix=3)
i_mn2, N2_max = iamax (Hs2)
i_mn2, N2_max

(16, -7364.83921638747)

(16, -7364.83921638747)

abs (N_max - N2_max) / abs (N_max) * 100

0.692423954462578

0.692423954462578

print_indexed (zip (extract_from (ef, ix=2), extract_from (ef2, ix=2)))

  0: (1105.25731545310, 1117.75947805641)
  1: (718.718275076535, 730.433083966749)
  2: (633.549069179356, 642.712212281765)
  3: (591.836581624291, 599.659925833765)
  4: (563.210575255919, 569.939041097995)
  5: (532.078102761199, 537.877263858203)
  6: (499.833997417863, 504.871945297586)
  7: (457.462570773090, 462.010220242996)
  8: (460.049569028232, 465.936679091990)
  9: (396.446344558523, 400.629917342073)
 10: (348.097591827964, 351.282570307938)
 11: (290.356874002145, 292.589544118509)
 12: (228.424684578299, 229.863494658540)
 13: (160.761022106667, 161.558055334346)
 14: (90.8772462296946, 91.2116793924809)
 15: (3.55068351909369, 3.60781468671729)
 16: (1154.17579305683, 1159.58861750160)
 17: (853.695300583230, 862.593437418486)
 18: (747.795880243230, 756.646484169387)
 19: (711.118170869699, 718.884903640265)
 20: (681.097023512360, 687.920073100633)
 21: (651.060726225303, 656.998046132936)
 22: (615.394614960553, 620.531364216977)
 23: (586.965033587022, 591.812465318944)
 24: (531.708795044128, 537.042437098731)
 25: (504.494571275030, 508.792050634829)
 26: (449.154007453208, 452.341242861567)
 27: (392.786148321149, 395.069943533253)
 28: (330.431049823004, 331.919009974338)
 29: (263.678632955008, 264.523152202011)
 30: (189.296188144868, 189.659068084037)
 31: (125.355397002250, 125.445475252240)
 32: (-921.333804637156, -932.858040817396)
 33: (-1156.83406554658, -1172.95050733064)
 34: (-1193.41502448136, -1209.96877878201)
 35: (-1172.70110956581, -1188.19287832588)
 36: (-1130.77383687013, -1144.70812671854)
 37: (-1076.92076702333, -1089.21693184347)
 38: (-1013.26580312721, -1024.16402508931)
 39: (-956.528823882090, -967.111508895205)
 40: (-870.299611246859, -880.423321345340)
 41: (-769.752752892900, -777.726474788654)
 42: (-663.493180699778, -669.431322348229)
 43: (-546.621456971245, -550.746389953271)
 44: (-419.513100841093, -422.114678597494)
 45: (-283.532941119409, -284.938027439446)
 46: (-132.741652063528, -133.318693392026)
 47: (-18.1640685269318, -18.3046386625432)

  0: (1105.25731545310, 1117.75947805641)
  1: (718.718275076535, 730.433083966749)
  2: (633.549069179356, 642.712212281765)
  3: (591.836581624291, 599.659925833765)
  4: (563.210575255919, 569.939041097995)
  5: (532.078102761199, 537.877263858203)
  6: (499.833997417863, 504.871945297586)
  7: (457.462570773090, 462.010220242996)
  8: (460.049569028232, 465.936679091990)
  9: (396.446344558523, 400.629917342073)
 10: (348.097591827964, 351.282570307938)
 11: (290.356874002145, 292.589544118509)
 12: (228.424684578299, 229.863494658540)
 13: (160.761022106667, 161.558055334346)
 14: (90.8772462296946, 91.2116793924809)
 15: (3.55068351909369, 3.60781468671729)
 16: (1154.17579305683, 1159.58861750160)
 17: (853.695300583230, 862.593437418486)
 18: (747.795880243230, 756.646484169387)
 19: (711.118170869699, 718.884903640265)
 20: (681.097023512360, 687.920073100633)
 21: (651.060726225303, 656.998046132936)
 22: (615.394614960553, 620.531364216977)
 23: (586.965033587022, 591.812465318944)
 24: (531.708795044128, 537.042437098731)
 25: (504.494571275030, 508.792050634829)
 26: (449.154007453208, 452.341242861567)
 27: (392.786148321149, 395.069943533253)
 28: (330.431049823004, 331.919009974338)
 29: (263.678632955008, 264.523152202011)
 30: (189.296188144868, 189.659068084037)
 31: (125.355397002250, 125.445475252240)
 32: (-921.333804637156, -932.858040817396)
 33: (-1156.83406554658, -1172.95050733064)
 34: (-1193.41502448136, -1209.96877878201)
 35: (-1172.70110956581, -1188.19287832588)
 36: (-1130.77383687013, -1144.70812671854)
 37: (-1076.92076702333, -1089.21693184347)
 38: (-1013.26580312721, -1024.16402508931)
 39: (-956.528823882090, -967.111508895205)
 40: (-870.299611246859, -880.423321345340)
 41: (-769.752752892900, -777.726474788654)
 42: (-663.493180699778, -669.431322348229)
 43: (-546.621456971245, -550.746389953271)
 44: (-419.513100841093, -422.114678597494)
 45: (-283.532941119409, -284.938027439446)
 46: (-132.741652063528, -133.318693392026)
 47: (-18.1640685269318, -18.3046386625432)

print_indexed (zip (extract_from (ef, ix=5), extract_from (ef2, ix=5)))

  0: (202.615529560617, 202.424956850654)
  1: (523.284996367237, 530.238295048884)
  2: (601.578442857064, 610.308852948246)
  3: (609.490534309890, 618.253837227884)
  4: (598.695734108930, 606.830862860334)
  5: (577.086769605483, 584.344986545879)
  6: (555.803232354139, 562.153804846334)
  7: (496.479254853853, 501.174829803219)
  8: (473.853266688340, 479.793404003283)
  9: (421.655161064949, 426.443904480724)
 10: (373.136306697634, 376.841778229721)
 11: (318.196772392925, 320.882895294738)
 12: (258.752078734423, 260.556623263143)
 13: (192.655694889701, 193.726348046967)
 14: (129.190968544439, 129.710878705311)
 15: (18.1640685269433, 18.3046386625382)
 16: (300.451363628612, 303.181638041868)
 17: (646.488929681218, 653.587786195461)
 18: (718.639109448411, 727.285519679074)
 19: (728.179714969467, 736.714817847970)
 20: (716.371669049083, 724.260146272026)
 21: (697.586903376460, 704.648444266027)
 22: (664.905657322525, 670.914556150274)
 23: (652.843143003277, 657.919052237408)
 24: (565.638371309532, 570.893286089834)
 25: (525.841425133200, 530.324114872722)
 26: (474.924596033293, 478.363610850815)
 27: (420.535207237905, 423.049602147708)
 28: (360.517188336770, 362.214658968308)
 29: (296.967151184179, 297.987911100151)
 30: (220.323097857868, 220.817894326387)
 31: (177.929851303346, 178.096091522082)
 32: (-1154.14666421184, -1165.77507546036)
 33: (-1394.28480992445, -1410.23427036485)
 34: (-1429.75728031811, -1446.17042331934)
 35: (-1409.27673848183, -1424.63489094859)
 36: (-1367.43239527438, -1381.25819240494)
 37: (-1312.98151833700, -1325.17980848302)
 38: (-1251.87069090954, -1262.72702146921)
 39: (-1184.55193804741, -1194.96148933613)
 40: (-1070.13294258456, -1079.68533672467)
 41: (-974.995432586404, -982.665357734280)
 42: (-867.710744354456, -873.433554384065)
 43: (-750.966257060916, -754.968612122050)
 44: (-624.195821291778, -626.737811170317)
 45: (-486.263339329045, -487.646979184181)
 46: (-345.678494860118, -346.263369578626)
 47: (-177.929851303351, -178.096091522086)

  0: (202.615529560617, 202.424956850654)
  1: (523.284996367237, 530.238295048884)
  2: (601.578442857064, 610.308852948246)
  3: (609.490534309890, 618.253837227884)
  4: (598.695734108930, 606.830862860334)
  5: (577.086769605483, 584.344986545879)
  6: (555.803232354139, 562.153804846334)
  7: (496.479254853853, 501.174829803219)
  8: (473.853266688340, 479.793404003283)
  9: (421.655161064949, 426.443904480724)
 10: (373.136306697634, 376.841778229721)
 11: (318.196772392925, 320.882895294738)
 12: (258.752078734423, 260.556623263143)
 13: (192.655694889701, 193.726348046967)
 14: (129.190968544439, 129.710878705311)
 15: (18.1640685269433, 18.3046386625382)
 16: (300.451363628612, 303.181638041868)
 17: (646.488929681218, 653.587786195461)
 18: (718.639109448411, 727.285519679074)
 19: (728.179714969467, 736.714817847970)
 20: (716.371669049083, 724.260146272026)
 21: (697.586903376460, 704.648444266027)
 22: (664.905657322525, 670.914556150274)
 23: (652.843143003277, 657.919052237408)
 24: (565.638371309532, 570.893286089834)
 25: (525.841425133200, 530.324114872722)
 26: (474.924596033293, 478.363610850815)
 27: (420.535207237905, 423.049602147708)
 28: (360.517188336770, 362.214658968308)
 29: (296.967151184179, 297.987911100151)
 30: (220.323097857868, 220.817894326387)
 31: (177.929851303346, 178.096091522082)
 32: (-1154.14666421184, -1165.77507546036)
 33: (-1394.28480992445, -1410.23427036485)
 34: (-1429.75728031811, -1446.17042331934)
 35: (-1409.27673848183, -1424.63489094859)
 36: (-1367.43239527438, -1381.25819240494)
 37: (-1312.98151833700, -1325.17980848302)
 38: (-1251.87069090954, -1262.72702146921)
 39: (-1184.55193804741, -1194.96148933613)
 40: (-1070.13294258456, -1079.68533672467)
 41: (-974.995432586404, -982.665357734280)
 42: (-867.710744354456, -873.433554384065)
 43: (-750.966257060916, -754.968612122050)
 44: (-624.195821291778, -626.737811170317)
 45: (-486.263339329045, -487.646979184181)
 46: (-345.678494860118, -346.263369578626)
 47: (-177.929851303351, -178.096091522086)

Mij = extract_from (ef, ix=2)
M2ij = extract_from (ef2, ix=2)
ij_mm, Mij_max = iamax (Mij)
ij_mm2, M2ij_max = iamax (M2ij)
print (ij_mm, Mij_max)
print (ij_mm2, M2ij_max)

(34, -1193.41502448136)
(34, -1209.96877878201)

(34, -1193.41502448136)
(34, -1209.96877878201)

Mji = extract_from (ef, ix=5)
M2ji = extract_from (ef2, ix=5)
ji_mm, Mji_max = iamax (Mji)
ji_mm2, M2ji_max = iamax (M2ji)
print (ji_mm, Mji_max)
print (ji_mm2, M2ji_max)

(34, -1429.75728031811)
(34, -1446.17042331934)

(34, -1429.75728031811)
(34, -1446.17042331934)

M2ji_max - Mji_max

-16.4131430012351

-16.4131430012351

(abs (M2ji_max - Mji_max)) / abs (Mji_max) * 100

1.14796708694383

1.14796708694383

NS_P-delta_(v0.8)

5 days ago by fresl