II. Ensamblaje y Solución (Algoritmo Unificado) 🧩
#

Una vez que tenemos la matriz de rigidez de cada elemento en ejes globales ($\bar{\mathbf{k}}^{(e)}$) para cada componente de la estructura (como la barra axial $4 \times 4$ vista en la Sección I), el siguiente paso es combinarlas para formar la matriz del sistema completo ($\mathbf{K}$) y luego resolverlo.

1. Ensamblaje de la Matriz Global ($\mathbf{K}$) por Superposición Directa
#

El corazón del Método de la Rigidez reside en el proceso de ensamblaje. Este método se basa en los principios de equilibrio nodal y compatibilidad de desplazamientos.

Principio: La Matriz de Rigidez Global $\mathbf{K}$ del sistema completo se obtiene sumando las contribuciones de rigidez de cada matriz $\bar{\mathbf{k}}^{(e)}$ individual en las posiciones correctas, según cómo los elementos se conectan en los nudos.
Grados de Libertad Globales (GDL): Antes de ensamblar, debemos numerar todos los posibles movimientos (desplazamientos $u_x, u_y$ y rotaciones $\theta_z$ si aplica) en cada nudo de la estructura. Este conjunto forma los GDL globales del sistema. Si hay $N_{nudos}$ y cada nudo tiene $g$ GDL (e.g., $g=2$ para armaduras 2D), la matriz $\mathbf{K}$ tendrá dimensiones $(N_{nudos} \times g) \times (N_{nudos} \times g)$. Por ejemplo, la estructura siguiente tiene 3 elementos y 3 nudos. Cada nudo tiene 2 GDL (Movimiento $u_x, u_y$), por lo que, la matriz $\mathbf{K}$ tendrá dimensiones $(3 \times 2) \times (3 \times 2) = (6 \times 6)$.

Convención de numeración: En el texto y en las figuras los GDL se numeran desde 1, como en la mayoría de los libros. En Python los índices parten en 0, por lo que el GDL $n$ del texto corresponde al índice n - 1 del código:
Nudo GDL (texto) Índices (código)
1 1 $(u_1)$, 2 $(v_1)$ 0, 1
2 3 $(u_2)$, 4 $(v_2)$ 2, 3
3 5 $(u_3)$, 6 $(v_3)$ 4, 5
En general, si cada nudo tiene $g$ GDL, los índices del nudo $i$ (numerado desde 1) en el código van desde g*(i-1) hasta g*(i-1) + g - 1.
Mapeo y Suma (Superposición):
1. Inicializar $\mathbf{K}$: Crear una matriz global $\mathbf{K}$ llena de ceros, con tamaño igual al número total de GDL globales.
2. Iterar por Elementos: Para cada elemento de la estructura:
  - Calcular su $\bar{\mathbf{k}}^{(e)}$ ($4 \times 4$ para barra axial 2D).
  - Identificar los índices de GDL globales asociados a los nudos de ese elemento. Por ejemplo, si una barra conecta el Nudo 1 (GDL globales 1, 2) y el Nudo 2 (GDL globales 3, 4).
  - Sumar: Añadir cada término $(i, j)$ de la $\bar{\mathbf{k}}^{(e)}$ del elemento al término correspondiente $(I, J)$ de la matriz $\mathbf{K}$ global, donde $I$ y $J$ son los índices globales mapeados desde los índices locales $i$ y $j$.
Ejemplo Conceptual: El proceso se puede desarrollar de la siguiente manera para la figura mostrada más arriba. En la figura, $K^1, K^2, K^3$ identifican a los elementos 1, 2 y 3. Para cada elemento $e$, $\mathbf{K}^{(e)}$ es su matriz $\bar{\mathbf{k}}^{(e)}$ expandida al tamaño de la estructura (con ceros en los GDL que no le corresponden), de modo que $\mathbf{K} = \mathbf{K}^{(1)} + \mathbf{K}^{(2)} + \mathbf{K}^{(3)}$.
1. Inicializar $\mathbf{K}$: $$ \mathbf{K} = \begin{array}{cc}% & \begin{matrix} u_1\;\; & v_1\;\; & u_2\;\; & v_2\;\; & u_3\;\; & v_3\;\; \end{matrix} \\ \begin{matrix} F_{u_1} \\ F_{v_1} \\ F_{u_2} \\ F_{v_2} \\ F_{u_3} \\ F_{v_3} \end{matrix} & \begin{pmatrix} K_{11} & K_{12} & K_{13} & K_{14} & K_{15} & K_{16}\\ K_{21} & K_{22} & K_{23} & K_{24} & K_{25} & K_{26}\\ K_{31} & K_{32} & K_{33} & K_{34} & K_{35} & K_{36}\\ K_{41} & K_{42} & K_{43} & K_{44} & K_{45} & K_{46}\\ K_{51} & K_{52} & K_{53} & K_{54} & K_{55} & K_{56}\\ K_{61} & K_{62} & K_{63} & K_{64} & K_{65} & K_{66} \end{pmatrix} % \end{array} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}$$
2.1. Iterar por Elementos (elemento 1):
$$\bar{\mathbf{k}}^{(1)} = \begin{array}{cc}% & \begin{matrix} u_1\;\; & v_1\;\; & u_2\;\; & v_2\;\; \end{matrix} \\ \begin{matrix} F_{u_1} \\ F_{v_1} \\ F_{u_2} \\ F_{v_2} \end{matrix} & \begin{pmatrix} \bar{k}_{11}^{(1)} & \bar{k}_{12}^{(1)} & \bar{k}_{13}^{(1)} & \bar{k}_{14}^{(1)}\\ \bar{k}_{21}^{(1)} & \bar{k}_{22}^{(1)} & \bar{k}_{23}^{(1)} & \bar{k}_{24}^{(1)}\\ \bar{k}_{31}^{(1)} & \bar{k}_{32}^{(1)} & \bar{k}_{33}^{(1)} & \bar{k}_{34}^{(1)}\\ \bar{k}_{41}^{(1)} & \bar{k}_{42}^{(1)} & \bar{k}_{43}^{(1)} & \bar{k}_{44}^{(1)}\\ \end{pmatrix} % \end{array} \implies$$$$ \mathbf{K}^{(1)} = \begin{array}{cc}% & \begin{matrix} u_1\;\; & v_1\;\; & u_2\;\; & v_2\;\; & u_3\;\; & v_3\;\; \end{matrix} \\ \begin{matrix} F_{u_1} \\ F_{v_1} \\ F_{u_2} \\ F_{v_2} \\ F_{u_3} \\ F_{v_3} \end{matrix} & \begin{pmatrix} \bar{k}_{11}^{(1)} & \bar{k}_{12}^{(1)} & \bar{k}_{13}^{(1)} & \bar{k}_{14}^{(1)} & 0 & 0\\ \bar{k}_{21}^{(1)} & \bar{k}_{22}^{(1)} & \bar{k}_{23}^{(1)} & \bar{k}_{24}^{(1)} & 0 & 0\\ \bar{k}_{31}^{(1)} & \bar{k}_{32}^{(1)} & \bar{k}_{33}^{(1)} & \bar{k}_{34}^{(1)} & 0 & 0\\ \bar{k}_{41}^{(1)} & \bar{k}_{42}^{(1)} & \bar{k}_{43}^{(1)} & \bar{k}_{44}^{(1)} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} % \end{array}$$
Se le suma $\mathbf{K}^{(1)}$ a la matriz $\mathbf{K}$
2.2. Iterar por Elementos (elemento 2):
$$\bar{\mathbf{k}}^{(2)} = \begin{array}{cc}% & \begin{matrix} u_2\;\; & v_2\;\; & u_3\;\; & v_3\;\; \end{matrix} \\ \begin{matrix} F_{u_2} \\ F_{v_2} \\ F_{u_3} \\ F_{v_3} \end{matrix} & \begin{pmatrix} \bar{k}_{11}^{(2)} & \bar{k}_{12}^{(2)} & \bar{k}_{13}^{(2)} & \bar{k}_{14}^{(2)}\\ \bar{k}_{21}^{(2)} & \bar{k}_{22}^{(2)} & \bar{k}_{23}^{(2)} & \bar{k}_{24}^{(2)}\\ \bar{k}_{31}^{(2)} & \bar{k}_{32}^{(2)} & \bar{k}_{33}^{(2)} & \bar{k}_{34}^{(2)}\\ \bar{k}_{41}^{(2)} & \bar{k}_{42}^{(2)} & \bar{k}_{43}^{(2)} & \bar{k}_{44}^{(2)}\\ \end{pmatrix} % \end{array} \implies$$$$ \mathbf{K}^{(2)} = \begin{array}{cc}% & \begin{matrix} u_1\;\; & v_1\;\; & u_2\;\; & v_2\;\; & u_3\;\; & v_3\;\; \end{matrix} \\ \begin{matrix} F_{u_1} \\ F_{v_1} \\ F_{u_2} \\ F_{v_2} \\ F_{u_3} \\ F_{v_3} \end{matrix} & \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 &\bar{k}_{11}^{(2)} & \bar{k}_{12}^{(2)} & \bar{k}_{13}^{(2)} & \bar{k}_{14}^{(2)} \\ 0 & 0 &\bar{k}_{21}^{(2)} & \bar{k}_{22}^{(2)} & \bar{k}_{23}^{(2)} & \bar{k}_{24}^{(2)} \\ 0 & 0 &\bar{k}_{31}^{(2)} & \bar{k}_{32}^{(2)} & \bar{k}_{33}^{(2)} & \bar{k}_{34}^{(2)} \\ 0 & 0 &\bar{k}_{41}^{(2)} & \bar{k}_{42}^{(2)} & \bar{k}_{43}^{(2)} & \bar{k}_{44}^{(2)} \\ \end{pmatrix} % \end{array}$$
Se le suma $\mathbf{K}^{(2)}$ a la matriz $\mathbf{K}$
2.3. Iterar por Elementos (elemento 3):
$$\bar{\mathbf{k}}^{(3)} = \begin{array}{cc}% & \begin{matrix} u_1\;\; & v_1\;\; & u_3\;\; & v_3\;\; \end{matrix} \\ \begin{matrix} F_{u_1} \\ F_{v_1} \\ F_{u_3} \\ F_{v_3} \end{matrix} & \begin{pmatrix} \bar{k}_{11}^{(3)} & \bar{k}_{12}^{(3)} & \bar{k}_{13}^{(3)} & \bar{k}_{14}^{(3)}\\ \bar{k}_{21}^{(3)} & \bar{k}_{22}^{(3)} & \bar{k}_{23}^{(3)} & \bar{k}_{24}^{(3)}\\ \bar{k}_{31}^{(3)} & \bar{k}_{32}^{(3)} & \bar{k}_{33}^{(3)} & \bar{k}_{34}^{(3)}\\ \bar{k}_{41}^{(3)} & \bar{k}_{42}^{(3)} & \bar{k}_{43}^{(3)} & \bar{k}_{44}^{(3)}\\ \end{pmatrix} % \end{array} \implies$$$$ \mathbf{K}^{(3)} = \begin{array}{cc}% & \begin{matrix} u_1\;\; & v_1\;\; & u_2\;\; & v_2\;\; & u_3\;\; & v_3\;\; \end{matrix} \\ \begin{matrix} F_{u_1} \\ F_{v_1} \\ F_{u_2} \\ F_{v_2} \\ F_{u_3} \\ F_{v_3} \end{matrix} & \begin{pmatrix} \bar{k}_{11}^{(3)} & \bar{k}_{12}^{(3)} & 0 & 0 & \bar{k}_{13}^{(3)} & \bar{k}_{14}^{(3)} \\ \bar{k}_{21}^{(3)} & \bar{k}_{22}^{(3)} & 0 & 0 & \bar{k}_{23}^{(3)} & \bar{k}_{24}^{(3)} \\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ \bar{k}_{31}^{(3)} & \bar{k}_{32}^{(3)} & 0 & 0 & \bar{k}_{33}^{(3)} & \bar{k}_{34}^{(3)} \\ \bar{k}_{41}^{(3)} & \bar{k}_{42}^{(3)} & 0 & 0 & \bar{k}_{43}^{(3)} & \bar{k}_{44}^{(3)} \\ \end{pmatrix} % \end{array}$$
Se le suma $\mathbf{K}^{(3)}$ a la matriz $\mathbf{K}$

Nudo	GDL (texto)	Índices (código)
1	1 \((u_1)\), 2 \((v_1)\)	0, 1
2	3 \((u_2)\), 4 \((v_2)\)	2, 3
3	5 \((u_3)\), 6 \((v_3)\)	4, 5

2. Aplicación de Condiciones de Borde (Restricciones) 🚫
#

La matriz $\mathbf{K}$ ensamblada representa la estructura antes de aplicar los soportes (apoyos). Para obtener una solución única y físicamente realista, debemos incorporar las condiciones de borde, que son los GDL con desplazamientos conocidos (generalmente cero en los apoyos).

Modelado de Apoyos:
- Apoyo Fijo (Pinned): $u_x = 0, u_y = 0$.
- Apoyo Móvil (Roller): Restringe el desplazamiento perpendicular a la superficie de rodadura (e.g., $u_y = 0$).
- Empotramiento (Fixed): $u_x = 0, u_y = 0, \theta_z = 0$ (para vigas/marcos).
Partición del Sistema: El sistema global $\mathbf{K} \mathbf{u} = \mathbf{f}$ se puede particionar conceptualmente separando los GDL libres (desplazamientos incógnitos $\mathbf{u}_L$) de los GDL restringidos (desplazamientos conocidos $\mathbf{u}_R$, usualmente $\mathbf{0}$).
$$ \begin{pmatrix} \mathbf{K}_{LL} & \mathbf{K}_{LR} \\ \mathbf{K}_{RL} & \mathbf{K}_{RR} \end{pmatrix} \begin{pmatrix} \mathbf{u}_L \\ \mathbf{u}_R \end{pmatrix} = \begin{pmatrix} \mathbf{f}_L \\ \mathbf{f}_R \end{pmatrix} $$
Donde:
- $\mathbf{u}_L$: Desplazamientos incógnitos (libres).
- $\mathbf{u}_R$: Desplazamientos conocidos (restringidos).
- $\mathbf{f}_L$: Fuerzas aplicadas en los GDL libres.
- $\mathbf{f}_R$: Fuerzas en los GDL restringidos, que incluyen las reacciones desconocidas.
Reducción de la Matriz: Para encontrar los desplazamientos incógnitos ($\mathbf{u}_L$), solo necesitamos la primera fila de ecuaciones de la matriz particionada:
$$\mathbf{K}_{LL} \mathbf{u}_L + \mathbf{K}_{LR} \mathbf{u}_R = \mathbf{f}_L$$
Como $\mathbf{u}_R$ es conocido (usualmente $\mathbf{u}_R = \mathbf{0}$), el sistema se reduce a:
$$\mathbf{K}_{LL} \mathbf{u}_L = \mathbf{f}_L - \mathbf{K}_{LR} \mathbf{u}_R$$
Si $\mathbf{u}_R = \mathbf{0}$, simplemente resolvemos:
$$\mathbf{K}_{LL} \mathbf{u}_L = \mathbf{f}_L$$
Llamamos a $\mathbf{K}_{LL}$ la matriz de rigidez reducida. La eliminación de filas y columnas correspondientes a los GDL restringidos es la implementación práctica de esta reducción.

3. Solución Numérica 💡
#

Una vez obtenido el sistema reducido $\mathbf{K}_{LL} \mathbf{u}_L = \mathbf{f}_L$, este es un sistema de ecuaciones lineales estándar de la forma $A x = b$.

Herramienta de Solución: La forma más eficiente y numéricamente estable de resolver este sistema es utilizando la función numpy.linalg.solve:

import numpy as np
# K_LL: Matriz de rigidez reducida (Numpy array)
# f_L: Vector de fuerzas en los GDL libres (Numpy array)

u_L = np.linalg.solve(K_LL, f_L)

Evitar la Inversa: Aunque matemáticamente $\mathbf{u}_L = \mathbf{K}_{LL}^{-1} \mathbf{f}_L$, calcular explícitamente la inversa ($\mathbf{K}_{LL}^{-1}$) es computacionalmente más costoso y menos estable numéricamente que usar métodos directos como la descomposición LU (implementada en solve). Por lo tanto, np.linalg.solve es la opción preferida.
Aprovechar la estructura de $\mathbf{K}_{LL}$: Si los apoyos impiden todo movimiento de cuerpo rígido, $\mathbf{K}_{LL}$ es simétrica y definida positiva. SciPy puede aprovecharlo usando la descomposición de Cholesky, que requiere cerca de la mitad de operaciones que LU:
```
from scipy.linalg import solve
u_L = solve(K_LL, f_L, assume_a="pos")
```
Si la estructura es inestable (un mecanismo, por ejemplo por faltar un apoyo), $\mathbf{K}_{LL}$ deja de ser definida positiva y esta llamada suele arrojar un error LinAlgError: es una forma útil de detectar apoyos mal definidos. En cambio, np.linalg.solve puede no avisar y devolver desplazamientos absurdos (por ejemplo, del orden de $10^{12}$ m), por lo que conviene revisar también np.linalg.cond(K_LL).

4. Ejercicio Demostrativo: Armadura Simple de 2 Barras 💻
#

Analicemos una armadura simple de 2 barras conectadas en un nudo y apoyadas en los otros extremos.

Las flechas de la figura indican el sentido positivo de las cargas. En el código se usa $F_y = -10\,000$ N, es decir, la carga vertical actúa hacia abajo.

Datos:

Nudo 1: $(0, 0)$, Apoyo Fijo ($U_{1x}=0, U_{1y}=0$). GDL 1, 2.
Nudo 2: $(L, H)$, Libre. Carga $(F_x, F_y)$. GDL 3, 4.
Nudo 3: $(2L, 0)$, Apoyo Fijo ($U_{3x}=0, U_{3y}=0$). GDL 5, 6.
Barra 1 (1-2): Longitud $L_1$, Ángulo $\theta_1$.
Barra 2 (3-2): Longitud $L_2$, Ángulo $\theta_2$.
Propiedades: Mismo $A, E$ para ambas barras.

Pasos:

Calcular $\bar{\mathbf{k}}^{(1)}$ para Barra 1 (GDL involucrados: 1, 2, 3, 4).
Calcular $\bar{\mathbf{k}}^{(2)}$ para Barra 2 (GDL involucrados: 5, 6, 3, 4).
Ensamblar $\mathbf{K}$ global (6x6): Sumar las contribuciones en las posiciones correspondientes. Notar que los GDL 3 y 4 reciben aportes de ambas barras.
Definir $\mathbf{f}$ global (6x1): Solo $f_3=F_x, f_4=F_y$ son no nulos.
Aplicar Condiciones de Borde: GDL restringidos son 1, 2, 5, 6 ($\mathbf{u}_R = \mathbf{0}$). GDL libres son 3, 4 ($\mathbf{u}_L = [U_{2x}, U_{2y}]^T$).
Reducir el sistema: Extraer $\mathbf{K}_{LL}$ ($2 \times 2$, correspondiente a GDL 3 y 4) y $\mathbf{f}_L$ ($2 \times 1$, correspondiente a GDL 3 y 4). $$\mathbf{K}_{LL} = \begin{pmatrix} K_{33} & K_{34} \\ K_{43} & K_{44} \end{pmatrix}$$ $$\mathbf{f}_L = \begin{pmatrix} F_x \\ F_y \end{pmatrix}$$
Resolver: Calcular $\mathbf{u}_L$ a partir de $\mathbf{K}_{LL} \mathbf{u}_L = \mathbf{f}_L$ usando np.linalg.solve.

import numpy as np

def calcular_k_elem(A, E, L, theta_grados):
    # (Función de la Sección I: matriz 4x4 del elemento en ejes globales)
    AE_L = (A * E) / L
    angulo_rad = np.radians(theta_grados)
    c = np.cos(angulo_rad)
    s = np.sin(angulo_rad)
    c2 = c**2
    s2 = s**2
    cs = c * s
    k_e = AE_L * np.array([
        [ c2,  cs, -c2, -cs],
        [ cs,  s2, -cs, -s2],
        [-c2, -cs,  c2,  cs],
        [-cs, -s2,  cs,  s2]
    ])
    return k_e

# --- Datos del Problema ---
L = 4.0 # m
H = 3.0 # m
A = 0.001 # m^2
E = 210e9 # Pa (Acero)
Fx = 5000 # N
Fy = -10000 # N

# Geometría y Propiedades de Elementos
# Barra 1 (Nudo 1 a Nudo 2)
L1 = np.sqrt(L**2 + H**2)
theta1_rad = np.arctan2(H, L)
theta1_deg = np.degrees(theta1_rad)
gdl_1 = [0, 1, 2, 3] # Índices en el código de los GDL 1, 2, 3, 4

# Barra 2 (Nudo 3 a Nudo 2)
L2 = np.sqrt(L**2 + H**2) 
theta2_rad = np.arctan2(H, -L) # Angulo desde Nudo 3 a Nudo 2
theta2_deg = np.degrees(theta2_rad)
gdl_2 = [4, 5, 2, 3] # Índices en el código de los GDL 5, 6, 3, 4

# --- Ensamblaje ---
num_gdl_total = 6 
K = np.zeros((num_gdl_total, num_gdl_total))

# Calcular la matriz de cada barra en ejes globales
k_elem_1 = calcular_k_elem(A, E, L1, theta1_deg)
k_elem_2 = calcular_k_elem(A, E, L2, theta2_deg)

# Ensamblar la barra 1
for i_local, i_global in enumerate(gdl_1):
    for j_local, j_global in enumerate(gdl_1):
        K[i_global, j_global] += k_elem_1[i_local, j_local]

# Ensamblar la barra 2
for i_local, i_global in enumerate(gdl_2):
    for j_local, j_global in enumerate(gdl_2):
        K[i_global, j_global] += k_elem_2[i_local, j_local]

# Forma compacta equivalente (np.ix_ selecciona el bloque de filas y columnas):
# K[np.ix_(gdl_1, gdl_1)] += k_elem_1
# K[np.ix_(gdl_2, gdl_2)] += k_elem_2

# --- Definir Vector de Fuerzas Global ---
f = np.zeros(num_gdl_total)
f[2] = Fx # Fuerza en GDL 3 (U_2x), índice 2
f[3] = Fy # Fuerza en GDL 4 (U_2y), índice 3

# --- Aplicar Condiciones de Borde y Reducir ---
# Índices libres: 2, 3 (GDL 3, 4: U_2x, U_2y)
# Índices restringidos: 0, 1, 4, 5 (GDL 1, 2, 5, 6)
gdl_libres = [2, 3]
gdl_restringidos = [0, 1, 4, 5]

# Extraer K_LL
K_LL = K[np.ix_(gdl_libres, gdl_libres)]

# Extraer f_L
f_L = f[gdl_libres]

# --- Resolver para Desplazamientos Incógnitos ---
u_L = np.linalg.solve(K_LL, f_L)

print("--- Resultados del Análisis ---")
print("Matriz de Rigidez Global K (Completa 6x6):")
# print(K) # Descomentar para verla completa
print("\nMatriz Reducida K_LL (2x2):")
print(K_LL)
print("\nVector de Fuerzas en GDL libres f_L (2x1):")
print(f_L)
print("\nDesplazamientos Incógnitos u_L (U_2x, U_2y):")
print(u_L) # Unidades: metros

Con tu Módulo `rigidez.py` 🧰
#

Agrega a tu rigidez.py las funciones de esta sección:

ensamblar hace la superposición de un elemento con np.ix_ (la forma compacta de los dos ciclos for).
ensamblar_armadura arma la matriz $\mathbf{K}$ completa a partir de las listas de nudos y barras.
resolver aplica las condiciones de borde y resuelve el sistema reducido. Acepta desplazamientos impuestos $\mathbf{u}_R$ (asentamientos) y, si $\mathbf{K}_{LL}$ es singular, avisa con un error en lugar de entregar resultados absurdos.

# ============================ Sección II ============================

def gdl_barra(i, j):
    """Índices de GDL (base 0) de una barra que une los nudos i y j (2 GDL por nudo)."""
    return [2*i, 2*i + 1, 2*j, 2*j + 1]


def ensamblar(K, k_elem, gdl):
    """Suma la matriz de un elemento en la matriz K de la estructura (modifica K)."""
    K[np.ix_(gdl, gdl)] += k_elem


def ensamblar_armadura(nudos, barras, A, E):
    """Matriz K de una armadura 2D.

    nudos: lista de coordenadas (x, y). barras: lista de pares (i, j) de nudos.
    A y E pueden ser un solo valor o una lista con un valor por barra.
    """
    A = np.broadcast_to(A, len(barras))
    E = np.broadcast_to(E, len(barras))
    K = np.zeros((2*len(nudos), 2*len(nudos)))
    for (i, j), A_e, E_e in zip(barras, A, E):
        ensamblar(K, k_barra(A_e, E_e, nudos[i], nudos[j]), gdl_barra(i, j))
    return K


def resolver(K, f, gdl_restringidos, u_R=None):
    """Resuelve K u = f imponiendo los desplazamientos de los GDL restringidos.

    u_R: desplazamientos impuestos en gdl_restringidos (por defecto, ceros).
    Retorna el vector u completo. Si K_LL es (casi) singular, la estructura es
    un mecanismo y se lanza un error en lugar de entregar resultados absurdos.
    """
    gdl_libres = [g for g in range(len(f)) if g not in gdl_restringidos]
    u = np.zeros(len(f))
    if u_R is not None:
        u[gdl_restringidos] = u_R
    K_LL = K[np.ix_(gdl_libres, gdl_libres)]
    K_LR = K[np.ix_(gdl_libres, gdl_restringidos)]
    if np.linalg.cond(K_LL) > 1e12:
        raise np.linalg.LinAlgError("K_LL es singular: revisa los apoyos (¿la estructura es un mecanismo?)")
    u[gdl_libres] = np.linalg.solve(K_LL, f[gdl_libres] - K_LR @ u[gdl_restringidos])
    return u

Con el módulo, el ejercicio demostrativo completo se reduce a unas pocas líneas:

import numpy as np
from rigidez import ensamblar_armadura, resolver

nudos = [(0, 0), (4, 3), (8, 0)]   # Nudos 1, 2 y 3 (m)
barras = [(0, 1), (2, 1)]          # Barra 1: nudos 1-2; Barra 2: nudos 3-2

K = ensamblar_armadura(nudos, barras, A=0.001, E=210e9)
f = np.zeros(6)
f[[2, 3]] = [5000, -10000]         # Fx y Fy en el Nudo 2 (GDL 3 y 4)

u = resolver(K, f, gdl_restringidos=[0, 1, 4, 5])
print("Desplazamientos del Nudo 2 (m):", u[[2, 3]])

5. Ejercicios Propuestos 📝
#

Concepto: Explica por qué la matriz $\mathbf{K}$ de una estructura grande es rala (sparse, con muchos ceros). ¿Qué términos de $\mathbf{K}$ son distintos de cero?
Mapeo: Para una armadura triangular con nudos 1 $(0,0)$, 2 $(L,0)$ y 3 $(L/2, H)$, numera los GDL globales (1 a 6) y escribe, para cada una de las 3 barras, el vector de índices que usarías en el código (base 0).
Cargas y apoyos: Para la armadura del ejercicio 2, con una fuerza vertical $-P$ en el nudo 3 y una fuerza horizontal $Q$ en el nudo 2, escribe el vector $\mathbf{f}$ (6x1). Si el nudo 1 es un apoyo fijo y el nudo 2 tiene un rodillo que restringe $U_{2y}$, ¿cuáles son los GDL libres y restringidos?
Ensamblaje manual: Dados $\bar{\mathbf{k}}^{(1)}$ ($4 \times 4$, GDL 1 a 4) y $\bar{\mathbf{k}}^{(2)}$ ($4 \times 4$, GDL 3 a 6), escribe cómo calcularías el término $K_{33}$ de la matriz $\mathbf{K}$ (6x6). ¿Con qué índices accederías a ese término en el código?
Reducción manual: Dada una matriz $\mathbf{K}$ (4x4) y $\mathbf{f}$ (4x1), si los GDL 1 y 3 están restringidos a cero (índices 0 y 2 en el código), escribe la $\mathbf{K}_{LL}$ (2x2) y el $\mathbf{f}_L$ (2x1) resultantes.
Programación: Explica por qué K[np.ix_(gdl, gdl)] += k_elem (función ensamblar de tu módulo) equivale a los dos ciclos for del ejercicio demostrativo. Luego usa ensamblar_armadura para armar la $\mathbf{K}$ de la armadura del ejercicio 2 (con $L = 4$ m y $H = 3$ m) y verifica que es simétrica.
Programación: Modifica el ejercicio demostrativo para que la Barra 2 tenga el doble de área ($2A$). ¿Cuál de los dos desplazamientos del Nudo 2 cambia más? ¿Por qué?
Programación (estructura inestable): Cambia el Nudo 3 del ejercicio demostrativo por un rodillo que solo restringe $U_{3y}$.
- Escribe los nuevos gdl_libres y gdl_restringidos y resuelve con np.linalg.solve. ¿Son razonables los desplazamientos obtenidos?
- Calcula np.linalg.cond(K_LL) y resuelve también con scipy.linalg.solve(K_LL, f_L, assume_a="pos"). ¿Qué ocurre en cada caso?
- Resuélvelo también con la función resolver de tu módulo. ¿Qué ocurre?
- ¿Qué significa físicamente que $\mathbf{K}_{LL}$ sea singular? Cuenta barras, reacciones y GDL para justificarlo.
Aplicación (asentamiento): Si el apoyo del Nudo 3 del ejercicio demostrativo sufre un asentamiento de $-0.005$ m en la dirección Y, $\mathbf{u}_R$ deja de ser cero. Usa $\mathbf{K}_{LL} \mathbf{u}_L = \mathbf{f}_L - \mathbf{K}_{LR} \mathbf{u}_R$ para calcular los nuevos desplazamientos del Nudo 2. Verifica tu resultado con resolver(K, f, gdl_restringidos, u_R=...).
Interpretación: Si el resultado u_L del ejercicio demostrativo fuera [0.001, -0.002] (en metros), ¿qué significaría físicamente para el Nudo 2? ¿Te parecen valores razonables para una armadura de acero de 8 m de luz con barras de 5 m?