Ensamblaje paso a paso¶

import numpy as np
from IPython.display import display, Markdown
import pandas as pd

def readin():
    nodes = np.loadtxt('files/sprnodes.txt', ndmin=2)
    mats = np.loadtxt('files/sprmater.txt', ndmin=2)
    elements = np.loadtxt('files/spreles.txt', ndmin=2, dtype=np.int)
    loads = np.loadtxt('files/sprloads.txt', ndmin=2)
    return nodes, mats, elements, loads

def eqcounter(nodes):
    nn = nodes.shape[0]
    IBC = np.zeros((nn, 1), dtype=np.integer)
    neq = 0
    for cont in range(nn):
        IBC[cont] = int(nodes[cont, 2])
        if IBC[cont] == 0:
            IBC[cont] = neq
            neq = neq + 1
    return neq, IBC

def DME(nodes, elements):
    nels = elements.shape[0]
    DME_mat = np.zeros((nels, 2), dtype=np.integer)
    neq, IBC = eqcounter(nodes)
    ndof = 2
    nnodes = 2
    for ele in range(nels):
        for node in range(nnodes):
            pos = elements[ele, node + 3]
            DME_mat[ele, node] = IBC[pos]
    return DME_mat, IBC, neq

def uelspring(kcof):
    """
    Esta rutina calcula la matriz de rigidez para
    un elemento tipo resorte conformado por 2 nodos.

    Parámetros
    ----------
    kcof : float
      Coeficiente de rigidez del resorte (>0).

    Retorna
    -------
    kloc : ndarray
      Matriz de coeficientes de rigidez local para
      el elemento tipo resorte (2, 2).


    """
    kloc = np.array([
        [kcof, -kcof],
        [-kcof, kcof]])
    return kloc

def assem_msg(el, kloc, kglob):
    """Muestra la matriz local y global para un elemento dado"""
    display(Markdown("### Elemento {:d}".format(el)))
    display(Markdown("Matriz local"))
    display(pd.DataFrame(kloc))
    display(Markdown("Matriz global luego de ensamblar"))
    display(pd.DataFrame(kglob))
    return None

def assembly(elements, mats, nodes, neq, DME_mat, steps=True):
    IELCON = np.zeros([2], dtype=np.integer)
    KG = np.zeros((neq, neq))
    nels = elements.shape[0]
    nnodes = 2
    ndof = 2
    for el in range(nels):
        elcoor = np.zeros((nnodes))
        im = np.int(elements[el, 2])
        par = mats[im]
        for j in range(nnodes):
            IELCON[j] = elements[el, j+3]
            elcoor[j] = nodes[IELCON[j], 1]
        kloc = uelspring(par[0])
        dme = DME_mat[el, :ndof]
        for row in range(ndof):
            glob_row = dme[row]
            if glob_row != -1:
                for col in range(ndof):
                    glob_col = dme[col]
                    if glob_col != -1:
                        KG[glob_row, glob_col] = KG[glob_row, glob_col] +\
                            kloc[row, col]
        if steps:
            assem_msg(el, kloc, KG)
    return KG

def loadasem(loads, IBC, neq, nl):
    """
    Ensambla el vector de cargas o de valores conocidos
    en el sistema global de ecuaciones.

    Parámetros
    ----------
    loads : ndarray
      Arreglo con las cargas impuestas a las partículas.
    IBC : ndarray (int)
      Arreglo que almacena el identificador de grado de libertad
      a cada partícula.
    neq : int
      Numero de ecuaciones en el sistema.
    nl : int
      Numero de cargas en el sistema.

    Retorna
    -------
    rhs_glob : ndarray
      Arreglo con las cargas impuestas a las partículas,
      rhs se refiere a lado derecho (del inglés right-hand-side).

    """
    rhs_glob = np.zeros((neq))
    for cont in range(nl):
        il = int(loads[cont, 0])
        ilx = IBC[il]
        if ilx != -1:
            rhs_glob[ilx] = loads[cont, 1]
    return rhs_glob

nodes, mats, elements, loads = readin()