"""
Feldbett Connector FEM Simulation - FreeCAD Python Script
==========================================================
Simuliert die Belastung von Connector 2 (Q+D+B) unter Last.

VORAUSSETZUNGEN:
  1. FreeCAD 0.20 oder neuer
  2. CalculiX installiert:
       Windows : meist mit FreeCAD mitgeliefert (prüfe Edit→Preferences→FEM→CalculiX)
       Ubuntu  : sudo apt install calculix-ccx
       Mac     : brew install calculix
  3. Das Connector-Script (feldbett_connectors_v2.py) muss VORHER
     ausgeführt worden sein — dieses Script baut auf dem Dokument
     "Feldbett_Connectors_v2" auf.
     ODER: Beide Scripts nacheinander in einer FreeCAD-Session ausführen.

ANWENDUNG:
  Makro → Makro ausführen → diese Datei wählen
  → FreeCAD öffnet automatisch die FEM-Ansicht
  → Solver starten: Doppelklick auf "SolverCcxTools" im Modellbaum
    → "Write .inp file" → "Run CalculiX"
  → Ergebnisse: Doppelklick auf "CCX_Results" → Pipeline aktivieren

ERGEBNISSE INTERPRETIEREN:
  Von-Mises-Spannung [N/mm²]:
    < Streckgrenze/Sicherheit  → grün, sicher
    > Streckgrenze             → rot, versagt
  Verformung [mm]:
    Zeigt wo der Connector nachgibt
"""

import FreeCAD as App
import FreeCADGui as Gui
import Part
import math
import os

# ============================================================
#  PARAMETER — müssen identisch zu feldbett_connectors_v2.py sein!
# ============================================================

L          = 35.0
BREITE_AA  = 70.0
MODULE     = 3

ROHR_DA    = 25.0
ROHR_WAND  =  2.0
HUELSE_DA  = 32.0
HUELSE_WAND=  3.0
EINSTECK_TIEFE = 45.0
SPIEL      =  0.4

LAST_KG    = 200.0
SICHERHEIT =   2.0
STRECKGRENZE = 235.0   # N/mm² (S235 Stahl)
E_MODUL    = 210000.0  # N/mm²
POISSON    = 0.3

# ============================================================
#  GEOMETRIE (identisch zu Connector-Script)
# ============================================================

hwA  = BREITE_AA / 2.0
hwQ  = L / 2.0
dZ   = hwA - hwQ
legH = math.sqrt(0.75 * L**2 - dZ**2)

def norm(v):
    l = math.sqrt(sum(x**2 for x in v))
    return tuple(x/l for x in v)

vQ  = (0.0, 0.0, 1.0)
vD  = norm(( L/2, -legH, -dZ))
vDr = norm((-L/2,  legH,  dZ))
vB  = norm((-L/2, -legH,  dZ))

ROHR_DI   = ROHR_DA - 2.0 * ROHR_WAND
HUELSE_DI = ROHR_DA + 2.0 * SPIEL

g = 9.81
F_gesamt        = LAST_KG * SICHERHEIT * g       # N, mit Sicherheit
F_pro_connector = F_gesamt / (MODULE * 4)         # N pro Connector 2
# An Connector 2 greifen D und B an — je 2 Stutzen
F_pro_stutzen   = F_pro_connector / 2.0           # N pro Stutzen

print("=" * 55)
print("FEM SIMULATION — CONNECTOR 2")
print("=" * 55)
print(f"  Auslegungslast    = {LAST_KG*SICHERHEIT:.0f} kg ({SICHERHEIT:.0f}× Sicherheit)")
print(f"  Kraft pro Stutzen = {F_pro_stutzen:.1f} N")
print(f"  Streckgrenze      = {STRECKGRENZE:.0f} N/mm²")
print(f"  Zul. Spannung     = {STRECKGRENZE/SICHERHEIT:.0f} N/mm²")
print()

# ============================================================
#  HILFSFUNKTIONEN
# ============================================================

def fv(t):
    return App.Vector(t[0], t[1], t[2])

def rotation_to_dir(direction):
    z = App.Vector(0, 0, 1)
    d = App.Vector(*direction).normalize()
    cross = z.cross(d)
    dot   = z.dot(d)
    if cross.Length > 1e-6:
        angle = math.degrees(math.acos(max(-1.0, min(1.0, dot))))
        return App.Rotation(cross, angle)
    elif dot < 0:
        return App.Rotation(App.Vector(1, 0, 0), 180)
    else:
        return App.Rotation()

def make_huelse_solid(direction, origin=(0,0,0)):
    """Solide Hülse (für FEM — keine Bohrung, vereinfacht für Vernetzung)."""
    outer = Part.makeCylinder(HUELSE_DA / 2.0, EINSTECK_TIEFE)
    inner = Part.makeCylinder(HUELSE_DI / 2.0, EINSTECK_TIEFE + 1.0)
    tube  = outer.cut(inner)
    rot = rotation_to_dir(direction)
    tube.Placement = App.Placement(App.Vector(*origin), rot)
    return tube

# ============================================================
#  GEOMETRIE FÜR FEM AUFBAUEN
#  Vereinfachter Connector 2: 4 Hülsen zusammengeschweißt
#  (FEM arbeitet besser mit einem zusammenhängenden Solid)
# ============================================================

print("Erzeuge Connector-Geometrie für FEM...")

O = (0, 0, 0)

h_Ql = make_huelse_solid((0, 0, -1), O)
h_Qr = make_huelse_solid((0, 0, +1), O)
h_Dr = make_huelse_solid(vDr, O)
h_Bv = make_huelse_solid(vB,  O)

# Alle zu einem Solid verschmelzen (wichtig für FEM — muss ein Body sein)
print("  Verschmelze Hülsen (kann einen Moment dauern)...")
try:
    connector_solid = h_Ql.fuse(h_Qr).fuse(h_Dr).fuse(h_Bv)
    # FEM braucht ein sauberes Solid ohne lose Faces
    connector_solid = connector_solid.removeSplitter()
    print("  ✓ Solid erstellt")
except Exception as e:
    print(f"  ⚠ Fuse-Fehler: {e}")
    print("    Verwende h_Ql als Fallback für Test")
    connector_solid = h_Ql

# ============================================================
#  FEM DOKUMENT AUFSETZEN
# ============================================================

doc_name = "Feldbett_FEM"
if doc_name in App.listDocuments():
    App.closeDocument(doc_name)
doc = App.newDocument(doc_name)

print("Erzeuge FEM-Struktur...")

# Geometrie-Objekt
geo_obj = doc.addObject("Part::Feature", "Connector2_Solid")
geo_obj.Shape = connector_solid

doc.recompute()

# ============================================================
#  FEM ANALYSIS
# ============================================================

try:
    import ObjectsFem
    import FemGui
except ImportError as e:
    print(f"FEHLER: FEM-Modul nicht verfügbar: {e}")
    print("Stelle sicher dass FreeCAD mit FEM-Workbench installiert ist.")
    raise

# Analysis-Container
analysis = ObjectsFem.makeAnalysis(doc, "Analysis")

# --- Material ---
mat_obj = ObjectsFem.makeMaterialSolid(doc, "Material_Stahl")
mat = mat_obj.Material
mat["Name"]                  = "Stahl_S235"
mat["YoungsModulus"]         = f"{E_MODUL:.0f} MPa"
mat["PoissonRatio"]          = f"{POISSON}"
mat["Density"]               = "7900 kg/m^3"
mat["UltimateTensileStrength"]= f"{STRECKGRENZE:.0f} MPa"
mat["YieldStrength"]         = f"{STRECKGRENZE:.0f} MPa"
mat_obj.Material = mat
mat_obj.References = [(geo_obj, "Solid1")]
analysis.addObject(mat_obj)

# --- Mesh (Vernetzung) ---
# Netfgen oder GMSH — FreeCAD nutzt intern Netgen
mesh_obj = doc.addObject("Fem::FemMeshShapeNetgenObject", "FEMMesh")
mesh_obj.Shape    = geo_obj
mesh_obj.MaxSize  = ROHR_DA / 2.0    # Elementgröße = halber Rohrdurchmesser
mesh_obj.MinSize  = ROHR_WAND        # Mindestgröße = Wandstärke
mesh_obj.Fineness = 3                # 1=sehr grob, 5=sehr fein (3=mittel, gut für Start)
mesh_obj.Optimize = True
mesh_obj.SecondOrder = True          # Quadratische Elemente = genauer
analysis.addObject(mesh_obj)

print(f"  Mesh: MaxSize={ROHR_DA/2:.1f}mm, Fineness=3 (mittel)")

# --- Einspannung (Fixed Constraint) ---
# Q-Stange an einem Ende festhalten (simuliert Befestigung am Bett-Rahmen)
# Face-Auswahl: Stirnfläche der Q-links-Hülse
fix = ObjectsFem.makeConstraintFixed(doc, "Einspannung_Q_Ende")
# Face muss manuell ausgewählt werden — wir setzen eine plausible Referenz
# (In der Praxis: in FreeCAD GUI die richtige Face anklicken)
fix.References = [(geo_obj, "Face1")]   # ← ggf. in GUI anpassen
analysis.addObject(fix)

# --- Kraft an D-Stutzen ---
force_D = ObjectsFem.makeConstraintForce(doc, "Last_D_Stutzen")
force_D.Force    = F_pro_stutzen        # N
force_D.Direction = (geo_obj, ["Edge1"]) # Richtung entlang D-Achse
force_D.Reversed  = False
# Kraftvektor in D-Richtung
force_D.DirectionVector = App.Vector(*vD)
force_D.References = [(geo_obj, "Face2")]  # ← Stirnfläche D-Stutzen, in GUI anpassen
analysis.addObject(force_D)

# --- Kraft an B-Stutzen ---
force_B = ObjectsFem.makeConstraintForce(doc, "Last_B_Stutzen")
force_B.Force    = F_pro_stutzen
force_B.DirectionVector = App.Vector(*vB)
force_B.References = [(geo_obj, "Face3")]  # ← Stirnfläche B-Stutzen, in GUI anpassen
analysis.addObject(force_B)

# --- Solver (CalculiX) ---
solver = ObjectsFem.makeSolverCalculixCcxTools(doc, "SolverCcxTools")
solver.AnalysisType  = "static"
solver.GeometricalNonlinearity = "linear"
solver.ThermoMechSteadyState   = False
solver.MatrixSolverType        = "default"
solver.IterationsControlParameterTimeUse = False
solver.SplitInputWriter = False

# CalculiX-Pfad automatisch suchen
ccx_paths = [
    "/usr/bin/ccx",                          # Linux apt
    "/usr/local/bin/ccx",                    # Linux/Mac brew
    "C:/Program Files/FreeCAD/bin/ccx.exe",  # Windows
    "C:/Program Files (x86)/FreeCAD/bin/ccx.exe",
]
for path in ccx_paths:
    if os.path.exists(path):
        solver.ccxBinaryPath = path
        print(f"  ✓ CalculiX gefunden: {path}")
        break
else:
    print("  ⚠ CalculiX nicht automatisch gefunden.")
    print("    → Manuell setzen: Edit → Preferences → FEM → CalculiX")

analysis.addObject(solver)

# ============================================================
#  ABSCHLUSS & ANLEITUNG
# ============================================================

doc.recompute()

# FEM-Workbench aktivieren
try:
    Gui.activateWorkbench("FemWorkbench")
    FemGui.setActiveAnalysis(analysis)
except Exception:
    pass

try:
    Gui.activeDocument().activeView().fitAll()
except Exception:
    pass

print()
print("=" * 55)
print("FEM-SETUP ABGESCHLOSSEN")
print("=" * 55)
print("""
NÄCHSTE SCHRITTE IN FREECAD:

1. FACES ANPASSEN (wichtig!):
   Die Einspannungen und Lasten sind auf "Face1/2/3"
   gesetzt — das sind Platzhalter.
   → Im Modellbaum: Doppelklick auf "Einspannung_Q_Ende"
     → Stirnfläche der Q-Hülse anklicken → OK
   → Gleiches für "Last_D_Stutzen" und "Last_B_Stutzen"

2. MESH ERZEUGEN:
   → Doppelklick auf "FEMMesh" im Modellbaum
   → "Mesh parameters" prüfen → OK
   → Vernetzung startet automatisch (dauert 10-30 Sek.)

3. SIMULATION STARTEN:
   → Doppelklick auf "SolverCcxTools"
   → Button "Write .inp file"
   → Button "Run CalculiX"
   → Warten (je nach Mesh 30 Sek. bis 5 Min.)

4. ERGEBNISSE ANZEIGEN:
   → Im Modellbaum erscheint "CCX_Results"
   → Menü: FEM → Postprocessing → Apply pipeline to result
   → Wähle "Von Mises Stress" oder "Displacement"
   → Farbskala zeigt Spannungsverteilung
""")
print(f"GRENZWERTE:")
print(f"  Zulässige Spannung : {STRECKGRENZE/SICHERHEIT:.0f} N/mm²  (grün im Plot)")
print(f"  Streckgrenze       : {STRECKGRENZE:.0f} N/mm²  (orange = kritisch)")
print(f"  Simulierte Last    : {LAST_KG*SICHERHEIT:.0f} kg  ({SICHERHEIT:.0f}× Sicherheit)")