# fem_vakwerk.py - het dakspant als stijfheidsmodel: K u = F (het idee van de eindige-elementenmethode)
import math
import numpy as np

h = math.tan(math.radians(30))
knopen = {"A": (0, 0), "D": (1.5, 1.5 * h), "C": (3, 3 * h), "E": (4.5, 1.5 * h), "B": (6, 0), "F": (2, 0), "G": (4, 0)}
staven = [("A", "D"), ("D", "C"), ("C", "E"), ("E", "B"), ("A", "F"), ("F", "G"), ("G", "B"),
          ("D", "F"), ("F", "C"), ("C", "G"), ("G", "E")]
E = 11e6                                                      # vuren C24, gemiddeld (kN/m2)
EA = {st: E * (0.045 * 0.095 if i < 7 else 0.045 * 0.070) for i, st in enumerate(staven)}   # randen en diagonalen (kN)
idx = {k: i for i, k in enumerate(knopen)}


def spantlasten(q_links, q_rechts, q_onder):
    return {"D": (0, -1.5 * q_links), "C": (0, -0.75 * (q_links + q_rechts)), "E": (0, -1.5 * q_rechts),
            "F": (0, -2.0 * q_onder), "G": (0, -2.0 * q_onder),
            "A": (0, -0.75 * q_links - q_onder), "B": (0, -0.75 * q_rechts - q_onder)}


def element(p, q):
    """Een staaf van p naar q: rek = b . (u_p, u_q) / L, dus stijfheid EA/L * b b^T met b = (-c, -s, c, s)."""
    d = np.subtract(knopen[q], knopen[p])
    L = np.linalg.norm(d)
    b = np.concatenate([-d / L, d / L])
    vg = [2 * idx[p], 2 * idx[p] + 1, 2 * idx[q], 2 * idx[q] + 1]   # de vier vrijheidsgraden van de staaf
    return L, b, vg


def los(vaste_vg, lasten):
    n = 2 * len(knopen)
    K, F = np.zeros((n, n)), np.zeros(n)
    for st in staven:                                         # 1. elementen maken en 2. assembleren
        L, b, vg = element(*st)
        K[np.ix_(vg, vg)] += EA[st] / L * np.outer(b, b)
    for knoop, f in lasten.items():
        F[2 * idx[knoop]:2 * idx[knoop] + 2] += f
    vrij = [i for i in range(n) if i not in vaste_vg]         # 3. randvoorwaarden: vaste vrijheidsgraden weg
    u = np.zeros(n)
    u[vrij] = np.linalg.solve(K[np.ix_(vrij, vrij)], F[vrij])  # 4. oplossen
    krachten = [EA[st] / element(*st)[0] * element(*st)[1] @ u[element(*st)[2]] for st in staven]   # 5. nabewerken
    return u, np.array(krachten), K @ u - F


q_b = 0.60 / math.cos(math.radians(30)) * 0.60                # eigengewicht dak (kN/m, horizontaal)
q_s, q_o = 0.40 * 0.60, 0.30 * 0.60                           # sneeuw en plafond
steunen = [0, 1, 2 * idx["B"] + 1]                            # A vast in x en y, B alleen in y

u, s1, R = los(steunen, spantlasten(1.35 * q_b + 1.5 * q_s, 1.35 * q_b + 1.5 * q_s, 1.35 * q_o))
print("rekenwaarden, twee steunen:  " + " ".join(f"{p}{q} {x:5.2f}" for (p, q), x in zip(staven, s1)))
u_k, _, _ = los(steunen, spantlasten(q_b + q_s, q_b + q_s, q_o))
print(f"doorbuiging (karakteristiek): nok {-1000 * u_k[2 * idx['C'] + 1]:.2f} mm, "
      f"onderrand in F {-1000 * u_k[2 * idx['F'] + 1]:.2f} mm, rol B schuift {1000 * u_k[2 * idx['B']]:.2f} mm")
u_w, s2, R2 = los(steunen + [2 * idx["F"] + 1], spantlasten(1.35 * q_b + 1.5 * q_s, 1.35 * q_b + 1.5 * q_s, 1.35 * q_o))
print("met tussenwand onder F:      " + " ".join(f"{p}{q} {x:5.2f}" for (p, q), x in zip(staven, s2)))
print(f"reacties met tussenwand: A {R2[1]:.2f}, F {R2[2 * idx['F'] + 1]:.2f}, B {R2[2 * idx['B'] + 1]:.2f} kN")