FEA-Driven Structural Optimization

import cadquery as cq
import numpy as np
from typing import Tuple

def parametric_bracket(
    base_width: float = 60.0,     # mm, base plate width (X)
    base_depth: float = 40.0,     # mm, base plate depth (Y)
    web_height: float = 50.0,     # mm, vertical web height (Z)
    wall_thickness: float = 4.0,  # mm, plate thickness
    fillet_radius: float = 4.0,   # mm, base-web junction fillet
    n_ribs: int = 2,              # number of stiffening ribs
    rib_thickness: float = None,  # mm, defaults to wall_thickness
    hole_diameter: float = 6.5,   # mm, mounting hole diameter (M6 clearance)
    hole_inset: float = 10.0,     # mm, hole center distance from edges
) -> cq.Workplane:
    """Generate a parametric L-bracket with ribs and fillets.

    The bracket sits at the origin with the base plate on the
    XY plane and the web extending upward along Z.

    Args:
        base_width: Width of the base plate (X direction)
        base_depth: Depth of the base plate (Y direction)
        web_height: Height of the vertical web (Z direction)
        wall_thickness: Thickness of plates
        fillet_radius: Fillet radius at base-web junction
        n_ribs: Number of triangular stiffening ribs
        rib_thickness: Rib thickness (defaults to wall_thickness)
        hole_diameter: Mounting hole diameter
        hole_inset: Distance from edge to hole center

    Returns:
        CadQuery solid of the bracket
    """
    t = wall_thickness
    if rib_thickness is None:
        rib_thickness = t

    # === Base plate ===
    base = (
        cq.Workplane("XY")
        .box(base_width, base_depth, t,
             centered=(True, True, False))
    )

    # === Vertical web ===
    # Web sits at the back edge of the base (Y = base_depth/2)
    web = (
        cq.Workplane("XY")
        .transformed(offset=(0, base_depth/2 - t/2, t))
        .box(base_width, t, web_height,
             centered=(True, True, False))
    )

    # Union base and web
    bracket = base.union(web)

    # === Fillet at base-web junction ===
    if fillet_radius > 0:
        # Select the interior edge where base meets web
        # Use NearestToPointSelector to target the concave junction edge
        try:
            bracket = (
                bracket
                .edges("|X")
                .edges(
                    cq.selectors.NearestToPointSelector(
                        (0, base_depth/2 - t, t)
                    )
                )
                .fillet(fillet_radius)
            )
        except Exception:
            pass  # Skip if fillet fails on complex geometry

    # === Stiffening ribs ===
    if n_ribs > 0:
        # Distribute ribs evenly across the width
        if n_ribs == 1:
            rib_positions = [0.0]
        else:
            margin = hole_inset + hole_diameter
            rib_span = base_width - 2 * margin
            rib_positions = np.linspace(
                -rib_span/2, rib_span/2, n_ribs
            ).tolist()

        rib_height = web_height * 0.7  # Ribs are 70% of web height
        rib_depth = base_depth * 0.6   # Ribs extend 60% of base depth

        for x_pos in rib_positions:
            # Triangular rib profile
            rib = (
                cq.Workplane("XY")
                .transformed(offset=(
                    x_pos, base_depth/2 - t, t
                ))
                .moveTo(0, 0)
                .lineTo(-rib_depth, 0)
                .lineTo(0, rib_height)
                .close()
                .extrude(rib_thickness)
                .translate((0, 0, 0))
            )

            # Center the rib extrusion on the x_pos
            rib = rib.translate((-rib_thickness/2, 0, 0))
            # Swap: extrude along X, so use workplane "YZ"
            rib = (
                cq.Workplane("YZ")
                .transformed(offset=(
                    x_pos - rib_thickness/2,
                    base_depth/2 - t,
                    t
                ))
                .moveTo(0, 0)
                .lineTo(-rib_depth, 0)
                .lineTo(0, rib_height)
                .close()
                .extrude(rib_thickness)
            )

            bracket = bracket.union(rib)

        # Fillet rib edges if radius permits
        if fillet_radius > 1:
            try:
                bracket = bracket.edges(
                    cq.selectors.BoxSelector(
                        (-base_width/2, -base_depth/2, t),
                        (base_width/2, base_depth/2, t + web_height)
                    )
                ).fillet(min(fillet_radius, rib_thickness * 0.4))
            except Exception:
                pass  # Skip if fillet fails on complex geometry

    # === Mounting holes in the base plate ===
    hole_positions = [
        (base_width/2 - hole_inset, base_depth/2 - hole_inset),
        (base_width/2 - hole_inset, -base_depth/2 + hole_inset),
        (-base_width/2 + hole_inset, base_depth/2 - hole_inset),
        (-base_width/2 + hole_inset, -base_depth/2 + hole_inset),
    ]

    for hx, hy in hole_positions:
        bracket = (
            bracket
            .faces("<Z")
            .workplane()
            .pushPoints([(hx, hy)])
            .hole(hole_diameter, t)
        )

    # === Load application hole in the web ===
    load_hole_z = t + web_height * 0.75  # 75% up the web
    bracket = (
        bracket
        .faces(">Y")
        .workplane()
        .pushPoints([(0, load_hole_z - t - web_height/2)])
        .hole(hole_diameter, t)
    )

    return bracket

# Generate a bracket with default parameters
bracket = parametric_bracket(
    wall_thickness=4.0,
    fillet_radius=4.0,
    n_ribs=2
)

# Export
cq.exporters.export(bracket, "bracket_t4_r2_f4.step")

import cadquery as cq

def bracket_mass(bracket: cq.Workplane,
                 density: float = 2.70e-6) -> float:
    """Calculate the mass of a bracket solid.

    Args:
        bracket: CadQuery solid
        density: Material density in kg/mm³
                 (Aluminum 6061: 2.70e-6 kg/mm³ = 2700 kg/m³)

    Returns:
        Mass in grams
    """
    # CadQuery volumes are in mm³
    volume_mm3 = bracket.val().Volume()
    mass_kg = volume_mm3 * density
    mass_g = mass_kg * 1000
    return mass_g

# Example
mass = bracket_mass(bracket)
print(f"Bracket volume: {bracket.val().Volume():.1f} mm³")
print(f"Bracket mass:   {mass:.2f} g")

import numpy as np
from itertools import product

def generate_parameter_space():
    """Define the 3D parameter sweep grid.

    Returns:
        List of (wall_thickness, n_ribs, fillet_radius) tuples
    """
    # Wall thickness: 2.0 to 6.0 mm in 1.0 mm steps
    thicknesses = [2.0, 3.0, 4.0, 5.0, 6.0]

    # Rib count: 0 to 4
    rib_counts = [0, 1, 2, 3, 4]

    # Fillet radius: 1.0 to 8.0 mm
    fillet_radii = [1.0, 2.0, 4.0, 6.0, 8.0]

    # Full factorial: 5 × 5 × 5 = 125 configurations
    space = list(product(thicknesses, rib_counts, fillet_radii))

    print(f"Parameter space: {len(space)} configurations")
    print(f"  Wall thickness: {thicknesses} mm")
    print(f"  Rib count:      {rib_counts}")
    print(f"  Fillet radius:  {fillet_radii} mm")

    return space

parameter_space = generate_parameter_space()

"""
FreeCAD FEM setup script.

Run this inside FreeCAD's Python console or via:
    freecad --console script.py

This script:
1. Imports a STEP file
2. Creates a FEM analysis
3. Assigns material properties
4. Creates mesh
5. Applies boundary conditions and loads
6. Solves with CalculiX
"""
import FreeCAD
import Part
import Fem
import ObjectsFem
import Mesh
import femmesh.femmesh2mesh

def setup_fem_analysis(
    step_file: str,
    output_dir: str,
    force_magnitude: float = 500.0,  # N
    mesh_size: float = 2.0,          # mm
) -> dict:
    """Set up and run a FEM analysis on an imported STEP bracket.

    Args:
        step_file: Path to the STEP file from CadQuery
        output_dir: Directory for solver output files
        force_magnitude: Applied force in Newtons
        mesh_size: Maximum mesh element size in mm

    Returns:
        Dictionary with analysis results
    """
    # --- Create new document ---
    doc = FreeCAD.newDocument("BracketFEM")

    # --- Import STEP geometry ---
    Part.insert(step_file, doc.Name)
    shape_obj = doc.Objects[0]
    doc.recompute()

    # --- Create FEM Analysis container ---
    analysis = ObjectsFem.makeAnalysis(doc, "Analysis")

    # --- Assign material: Aluminum 6061-T6 ---
    material = ObjectsFem.makeMaterialSolid(doc, "Aluminum6061T6")
    mat_dict = material.Material
    mat_dict["Name"] = "Aluminum 6061-T6"
    mat_dict["YoungsModulus"] = "68900 MPa"
    mat_dict["PoissonRatio"] = "0.33"
    mat_dict["Density"] = "2700 kg/m^3"
    mat_dict["ThermalExpansionCoefficient"] = "23.6 µm/m/K"
    mat_dict["YieldStrength"] = "276 MPa"
    mat_dict["UltimateTensileStrength"] = "310 MPa"
    material.Material = mat_dict
    analysis.addObject(material)

    # --- Create mesh with Gmsh ---
    femmesh = ObjectsFem.makeMeshGmsh(doc, "FEMMeshGmsh")
    femmesh.Part = shape_obj
    femmesh.CharacteristicLengthMax = mesh_size
    femmesh.CharacteristicLengthMin = mesh_size / 4
    femmesh.ElementOrder = "2nd"  # Quadratic elements
    femmesh.Algorithm2D = "Frontal"
    femmesh.Algorithm3D = "Delaunay"

    # Run Gmsh meshing
    from femmesh.gmshtools import GmshTools
    gmsh = GmshTools(femmesh)
    gmsh.create_mesh()

    analysis.addObject(femmesh)
    print(f"Mesh: {femmesh.FemMesh.NodeCount} nodes, "
          f"{femmesh.FemMesh.VolumeCount} elements")

    # --- Fixed constraint on base plate bottom face ---
    fixed = ObjectsFem.makeConstraintFixed(doc, "FixedBase")
    # Select the bottom face (Z = 0)
    bottom_faces = []
    for i, face in enumerate(shape_obj.Shape.Faces):
        com = face.CenterOfMass
        if abs(com.z) < 0.1:  # Bottom face
            bottom_faces.append(
                ("Face" + str(i + 1), face)
            )

    if bottom_faces:
        fixed.References = [
            (shape_obj, bottom_faces[0][0])
        ]
    analysis.addObject(fixed)

    # --- Force load on the web hole ---
    force = ObjectsFem.makeConstraintForce(doc, "AppliedForce")
    force.Force = force_magnitude  # Newtons

    # Select the cylindrical face of the load hole
    # (the hole near the top of the web)
    load_faces = []
    for i, face in enumerate(shape_obj.Shape.Faces):
        com = face.CenterOfMass
        # Load hole is at ~75% web height, on the far Y face
        if com.z > 30 and abs(com.y) > 15:
            if face.Surface.TypeId == "Part::GeomCylinder":
                load_faces.append(("Face" + str(i + 1), face))

    if load_faces:
        force.References = [
            (shape_obj, load_faces[0][0])
        ]
    force.Direction = (shape_obj, "Face1")  # Will be corrected
    force.Reversed = False
    analysis.addObject(force)

    # --- Solver: CalculiX ---
    solver = ObjectsFem.makeSolverCalculix(doc, "CalculiX")
    solver.WorkingDir = output_dir
    solver.AnalysisType = "static"
    solver.GeometricalNonlinearity = "linear"
    solver.MaterialNonlinearity = "linear"
    solver.EigenmodesCount = 0
    solver.IterationsMaximum = 2000
    analysis.addObject(solver)

    # --- Run solver ---
    doc.recompute()

    from femtools import ccxtools
    fea = ccxtools.FemToolsCcx(analysis, solver)
    fea.update_objects()
    fea.setup_working_dir()
    fea.setup_ccx()

    message = fea.check_prerequisites()
    if message:
        print(f"Prerequisites check: {message}")
        return {"error": message}

    fea.write_inp_file()
    fea.ccx_run()
    fea.load_results()

    doc.recompute()

    return doc, analysis, solver

# Example usage (inside FreeCAD):
# doc, analysis, solver = setup_fem_analysis(
#     "bracket_t4_r2_f4.step",
#     "/tmp/bracket_fem/",
#     force_magnitude=500.0,
#     mesh_size=2.0
# )

"""
Extract results from CalculiX .frd output file.

The .frd file format contains nodal displacements, stresses,
and strains computed by CalculiX.
"""
import os
import re
import numpy as np

def extract_frd_results(frd_path: str) -> dict:
    """Parse a CalculiX .frd file to extract key results.

    Args:
        frd_path: Path to the .frd results file

    Returns:
        Dictionary with max_von_mises, max_displacement, etc.
    """
    if not os.path.exists(frd_path):
        raise FileNotFoundError(f"FRD file not found: {frd_path}")

    von_mises_stresses = []
    displacements = []

    # Parse the .frd file
    with open(frd_path, 'r') as f:
        lines = f.readlines()

    # State machine for parsing
    in_stress_block = False
    in_disp_block = False

    for i, line in enumerate(lines):
        # Detect stress results block
        if "STRESS" in line and "100CL" in line:
            in_stress_block = True
            in_disp_block = False
            continue

        # Detect displacement block
        if "DISP" in line and "100CL" in line:
            in_disp_block = True
            in_stress_block = False
            continue

        # End of block
        if line.startswith(" -3"):
            in_stress_block = False
            in_disp_block = False
            continue

        # Parse stress data (nodal von Mises)
        if in_stress_block and line.startswith(" -1"):
            parts = line.split()
            if len(parts) >= 7:
                try:
                    # Columns: node, sxx, syy, szz, sxy, syz, sxz
                    sxx = float(parts[1])
                    syy = float(parts[2])
                    szz = float(parts[3])
                    sxy = float(parts[4])
                    syz = float(parts[5])
                    sxz = float(parts[6])

                    # Von Mises stress calculation
                    vm = np.sqrt(0.5 * (
                        (sxx - syy)**2 +
                        (syy - szz)**2 +
                        (szz - sxx)**2 +
                        6 * (sxy**2 + syz**2 + sxz**2)
                    ))
                    von_mises_stresses.append(vm)
                except (ValueError, IndexError):
                    continue

        # Parse displacement data
        if in_disp_block and line.startswith(" -1"):
            parts = line.split()
            if len(parts) >= 4:
                try:
                    dx = float(parts[1])
                    dy = float(parts[2])
                    dz = float(parts[3])
                    disp = np.sqrt(dx**2 + dy**2 + dz**2)
                    displacements.append(disp)
                except (ValueError, IndexError):
                    continue

    if not von_mises_stresses:
        return {"error": "No stress results found in .frd file"}

    results = {
        "max_von_mises_MPa": max(von_mises_stresses),
        "mean_von_mises_MPa": np.mean(von_mises_stresses),
        "max_displacement_mm": max(displacements) if displacements else 0,
        "n_stress_nodes": len(von_mises_stresses),
        "n_disp_nodes": len(displacements),
    }

    return results

# Alternative: Extract results from FreeCAD result object
def extract_freecad_results(doc) -> dict:
    """Extract results from FreeCAD FEM result objects.

    Args:
        doc: FreeCAD document with completed analysis

    Returns:
        Dictionary with stress and displacement results
    """
    # Find the result object
    result_obj = None
    for obj in doc.Objects:
        if hasattr(obj, "vonMises"):
            result_obj = obj
            break

    if result_obj is None:
        return {"error": "No result object found"}

    von_mises = np.array(result_obj.vonMises)
    disp_x = np.array(result_obj.DisplacementVectors)

    # Compute displacement magnitudes
    displacements = np.sqrt(
        disp_x[0::3]**2 + disp_x[1::3]**2 + disp_x[2::3]**2
    )

    return {
        "max_von_mises_MPa": float(np.max(von_mises)),
        "mean_von_mises_MPa": float(np.mean(von_mises)),
        "p95_von_mises_MPa": float(np.percentile(von_mises, 95)),
        "max_displacement_mm": float(np.max(displacements)),
        "n_nodes": len(von_mises),
    }

# Material: Aluminum 6061-T6
AL_6061_T6 = {
    "name": "Aluminum 6061-T6",
    "E_MPa": 68_900,           # Young's modulus
    "nu": 0.33,                # Poisson's ratio
    "rho_kg_m3": 2_700,        # Density
    "sigma_y_MPa": 276,        # Yield strength (0.2% offset)
    "sigma_ult_MPa": 310,      # Ultimate tensile strength
    "fatigue_limit_MPa": 96,   # Endurance limit at 5×10^8 cycles
    "alpha_per_C": 23.6e-6,    # Thermal expansion coefficient
}

# Design constraints for this optimization
DESIGN_CONSTRAINTS = {
    "max_von_mises_MPa": 276 / 2.0,  # σ_y / safety factor 2.0
    "max_displacement_mm": 0.5,        # Stiffness requirement
    "applied_force_N": 500,            # Load case
    "force_direction": (0, 0, -1),     # Downward (-Z)
}

print(f"Allowable stress: {DESIGN_CONSTRAINTS['max_von_mises_MPa']:.0f} MPa")
print(f"Maximum deflection: {DESIGN_CONSTRAINTS['max_displacement_mm']} mm")

import numpy as np
from typing import List, Dict, Tuple

def find_pareto_front(
    results: List[Dict],
    obj1_key: str = "mass_g",
    obj2_key: str = "max_von_mises_MPa",
    minimize_both: bool = True,
) -> List[int]:
    """Identify Pareto-optimal designs from a set of results.

    Args:
        results: List of result dictionaries
        obj1_key: Key for objective 1 (mass)
        obj2_key: Key for objective 2 (stress)
        minimize_both: Whether both objectives should be minimized

    Returns:
        List of indices into 'results' that are Pareto-optimal
    """
    # Filter to successful results only
    valid = [
        (i, r) for i, r in enumerate(results)
        if r.get("status") == "success"
    ]

    if not valid:
        return []

    n = len(valid)
    is_dominated = [False] * n

    for i in range(n):
        if is_dominated[i]:
            continue
        idx_i, ri = valid[i]
        m_i = ri[obj1_key]
        s_i = ri[obj2_key]

        for j in range(n):
            if i == j or is_dominated[j]:
                continue
            idx_j, rj = valid[j]
            m_j = rj[obj1_key]
            s_j = rj[obj2_key]

            # Does j dominate i?
            if minimize_both:
                if m_j <= m_i and s_j <= s_i and (
                    m_j < m_i or s_j < s_i
                ):
                    is_dominated[i] = True
                    break

    pareto_indices = [
        valid[i][0] for i in range(n) if not is_dominated[i]
    ]

    return pareto_indices

def select_best_design(
    results: List[Dict],
    pareto_indices: List[int],
    max_stress_MPa: float = 138.0,  # σ_y / 2.0
    max_displacement_mm: float = 0.5,
) -> int:
    """Select the best design from the Pareto front.

    Strategy: among Pareto-optimal designs that meet all
    constraints, select the one with minimum mass.

    Args:
        results: All results
        pareto_indices: Indices of Pareto-optimal designs
        max_stress_MPa: Maximum allowable stress
        max_displacement_mm: Maximum allowable displacement

    Returns:
        Index of the best design
    """
    feasible = []
    for idx in pareto_indices:
        r = results[idx]
        if r["max_von_mises_MPa"] <= max_stress_MPa:
            feasible.append(idx)

    if not feasible:
        print("WARNING: No Pareto-optimal design meets all "
              "constraints. Returning lowest-stress design.")
        return min(pareto_indices,
                   key=lambda i: results[i]["max_von_mises_MPa"])

    # Among feasible Pareto designs, pick lightest
    best = min(feasible, key=lambda i: results[i]["mass_g"])
    return best

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np

def plot_pareto_front(
    results: List[Dict],
    pareto_indices: List[int],
    best_index: int = None,
    constraint_stress: float = 138.0,
    title: str = "Bracket Optimization: Mass vs. Max Stress"
):
    """Plot the Pareto front with all design points.

    Args:
        results: All result dictionaries
        pareto_indices: Indices of Pareto-optimal designs
        best_index: Index of the selected best design
        constraint_stress: Maximum allowable stress line
    """
    fig, ax = plt.subplots(1, 1, figsize=(10, 7))

    # Separate valid results
    valid = [
        (i, r) for i, r in enumerate(results)
        if r.get("status") == "success"
    ]

    masses = [r["mass_g"] for _, r in valid]
    stresses = [r["max_von_mises_MPa"] for _, r in valid]
    indices = [i for i, _ in valid]

    # Color by rib count for insight
    rib_counts = [r["n_ribs"] for _, r in valid]
    scatter = ax.scatter(
        masses, stresses,
        c=rib_counts, cmap='viridis',
        s=40, alpha=0.5, edgecolors='gray',
        linewidth=0.5, zorder=2
    )
    cbar = plt.colorbar(scatter, ax=ax, label='Rib Count')

    # Highlight Pareto front
    pareto_masses = [results[i]["mass_g"] for i in pareto_indices]
    pareto_stresses = [
        results[i]["max_von_mises_MPa"] for i in pareto_indices
    ]

    # Sort Pareto points by mass for a clean line
    pareto_sorted = sorted(
        zip(pareto_masses, pareto_stresses)
    )
    pm_sorted = [p[0] for p in pareto_sorted]
    ps_sorted = [p[1] for p in pareto_sorted]

    ax.plot(pm_sorted, ps_sorted, 'r-o',
            markersize=8, linewidth=2,
            label='Pareto front', zorder=3)

    # Highlight best design
    if best_index is not None:
        bm = results[best_index]["mass_g"]
        bs = results[best_index]["max_von_mises_MPa"]
        ax.plot(bm, bs, 'r*', markersize=20, zorder=5)
        ax.annotate(
            f'BEST DESIGN\n'
            f't={results[best_index]["t"]}mm, '
            f'ribs={results[best_index]["n_ribs"]}, '
            f'fillet={results[best_index]["r_f"]}mm\n'
            f'mass={bm:.1f}g, σ={bs:.0f} MPa',
            xy=(bm, bs),
            xytext=(bm + 5, bs + 20),
            fontsize=9,
            arrowprops=dict(arrowstyle='->', color='red',
                            linewidth=2),
            bbox=dict(boxstyle='round,pad=0.4',
                      facecolor='lightyellow',
                      edgecolor='red', alpha=0.9),
            zorder=6
        )

    # Constraint line
    ax.axhline(y=constraint_stress, color='orange',
               linestyle='--', linewidth=1.5,
               label=f'Stress limit ({constraint_stress} MPa)')

    # Feasible region shading
    ax.axhspan(0, constraint_stress, alpha=0.05,
               color='green', label='Feasible region')

    # Labels and formatting
    ax.set_xlabel('Mass (g)', fontsize=13)
    ax.set_ylabel('Max von Mises Stress (MPa)', fontsize=13)
    ax.set_title(title, fontsize=14, fontweight='bold')
    ax.legend(loc='upper right', fontsize=10)
    ax.grid(True, alpha=0.3)

    # Add design direction arrows
    ax.annotate('', xy=(min(masses)*0.95, min(stresses)*0.9),
                xytext=(max(masses)*0.5, max(stresses)*0.5),
                arrowprops=dict(arrowstyle='->',
                                color='green', linewidth=2))
    ax.text(min(masses)*1.05, min(stresses)*0.95,
            'Ideal\ndirection',
            fontsize=9, color='green', style='italic')

    plt.tight_layout()
    plt.savefig("pareto_front.png", dpi=150, bbox_inches='tight')
    plt.show()

    return fig

# Usage:
# pareto_idx = find_pareto_front(all_results)
# best_idx = select_best_design(all_results, pareto_idx)
# plot_pareto_front(all_results, pareto_idx, best_idx)

import pandas as pd

def results_summary_table(
    results: List[Dict],
    pareto_indices: List[int],
    best_index: int
) -> pd.DataFrame:
    """Create a summary table of all valid results.

    Highlights Pareto-optimal and best designs.
    """
    rows = []
    for i, r in enumerate(results):
        if r.get("status") != "success":
            continue
        rows.append({
            "Config #": i,
            "t (mm)": r["t"],
            "Ribs": r["n_ribs"],
            "Fillet (mm)": r["r_f"],
            "Mass (g)": round(r["mass_g"], 1),
            "σ_max (MPa)": round(r["max_von_mises_MPa"], 1),
            "Pareto?": "Yes" if i in pareto_indices else "",
            "Best?": ">>> BEST <<<" if i == best_index else "",
        })

    df = pd.DataFrame(rows)

    # Sort by mass
    df = df.sort_values("Mass (g)").reset_index(drop=True)

    return df

# Usage:
# df = results_summary_table(all_results, pareto_idx, best_idx)
# print(df.to_string())
# df.to_csv("optimization_results.csv", index=False)