ch4_examples.md

Chapter 4: Worked Examples

---

Example 1: Curved Beam Optimisation

This example corresponds to examples/demo_beam.py. An S-shaped beam is imported from an IGES file, clamped at one end and guided at the other, then thickness-optimised.

Problem Setup

E         = 69e9        # Ti-6Al-4V, Young's modulus (Pa)
ultimate  = 900e6       # Ultimate tensile strength (Pa)
rho       = 4820        # Density (kg/m³)
beam_width = 0.008      # 8 mm out-of-plane width
initial_thickness = 0.003   # 3 mm starting thickness
Fx, Fy = 0, -333       # 333 N vertical load at free end

Geometry Import

The curve is imported from a SolidWorks-exported .igs file and resampled to 250 uniform-arc-length elements:

from feaopt import import_igs_curve

result = import_igs_curve(import_file, num_elements=250,
                          scale=0.001, flip=False, plot=True, verbose=True)

all_nodes     = result['nodes']
connectivity  = result['connectivity']
num_nodes     = result['num_nodes']
num_elements  = result['num_elements']

Boundary Conditions

Node 0 is fully clamped (all three DOFs fixed). The free end (last node) is constrained in x and θ but free in y (a "guided" support):

fixed_node = 0
free_node  = num_nodes - 1

constraints = np.array([
    [fixed_node, 0, 0],   # ux = 0
    [fixed_node, 1, 0],   # uy = 0
    [fixed_node, 2, 0],   # theta = 0
    [free_node,  0, 0],   # guided: ux = 0
    [free_node,  2, 0],   # guided: theta = 0
])

Analysis and Optimisation

A nonlinear analysis is run first to establish baseline deflection and stress. The optimiser then iterates with a safety factor of 2:

from feaopt import optimize_thickness

opt_options = {
    'alpha': 0.4,
    'safety_factor': 2,
    'max_iterations': 200,
    't_min': 0.0005,
    't_max': 0.020,
    'use_nonlinear': True,
    'verbose': True,
}

optimized_geom, opt_results = optimize_thickness(geometry, analysis, load_case, opt_options)

Typical Results

Metric	Initial	Optimised
Mass	7–10 g	3–5 g
Max displacement	0.5–2 mm	1–4 mm
Min factor of safety	≫ 2	≈ 2
Max stress	50–100 MPa	400–450 MPa
Thickness range	3 mm (uniform)	0.5–8 mm

The optimiser redistributes material from low-stress regions (near the neutral axis of curves) to high-stress regions (the outer fibres of bends), producing a variable-thickness organic profile.

---

Example 2: Truss Topology Optimisation

This example corresponds to examples/demo_truss.py. A Warren truss is defined programmatically, subdivided into fine elements, and optimised.

Problem Setup

bay_width  = 0.060      # 60 mm per bay
height     = 0.050      # 50 mm depth
n_bays     = 4          # 5 bottom nodes, 4 top nodes

E          = 200e9      # Steel
ultimate   = 400e6
rho        = 7850
beam_width = 0.008      # 8 mm width
initial_thickness = 0.003

M_tip  = 10             # 10 Nm tip moment
Fy_tip = -3000          # 3 kN downward tip force

num_sub = 50            # sub-elements per member

Truss Geometry

Bottom-chord nodes are spaced uniformly at y = 0. Top-chord nodes are at half-bay offsets at y = h:

bot_x = np.arange(n_bays + 1) * bay_width
warren_x = (np.arange(n_bays) + 0.5) * bay_width

Connectivity includes bottom chord, top chord, and diagonal members forming the characteristic alternating-V Warren pattern.

Subdivision

Each of the original members is split into 50 sub-elements:

from feaopt import subdivide_truss

sub_nodes, sub_conn, member_map = subdivide_truss(truss_nodes, conn, num_sub)

This creates approximately 843 nodes and 850 elements from the original 10 nodes and 17 members.

Optimisation

The optimiser uses linear analysis (sufficient for a stiff truss) with bounds that allow elements to thin:

opt_options = {
    'alpha': 0.35,
    'safety_factor': 2,
    'max_iterations': 300,
    't_min': 0.002,
    't_max': 0.02,
    'convergence_tol': 0.001,
    'move_limit': 0.15,
    'use_nonlinear': False,
    'verbose': True,
}

optimized_geom, opt_results = optimize_thickness(
    geometry, analysis, load_case, opt_options)

Interpreting Results

After optimisation, material concentrates along the primary load paths (top and bottom chords, and the most heavily loaded diagonals) while non-structural diagonals thin toward the minimum bound. The result resembles a continuous topology-optimised shape rather than a discrete truss.

The min_thickness_frac option in export_profile filters out vanished elements, leaving clean boundaries for CAD export:

from feaopt import export_profile

export_opts = {
    'min_thickness_frac': 0.15,
    'resolution': 5000,
    'export_filename': os.path.join(profile_dir, 'simple_truss'),
}
profile = export_profile(optimized_geom, export_opts)

Performance note: The default resolution of 50,000 pixels/m can produce very large images for truss geometries spanning hundreds of millimetres. Reducing it to 5,000 keeps the morphological operations fast while still producing smooth profiles.

---

Example 3: Minimal Custom Problem

A two-element L-bracket, built entirely from code with no file imports:

import numpy as np
from feaopt import StructureGeometry, FEAAnalysis, LoadCase

# Geometry
nodes = np.array([[0, 0], [0.1, 0], [0.1, 0.08]])
conn  = np.array([[0, 1], [1, 2]])

props = {
    'E': 200e9,
    'thickness': 0.005,
    'width': 0.012,
    'ultimate': 350e6,
    'rho': 7850,
}

constraints = np.array([[0, 0, 0], [0, 1, 0], [0, 2, 0]])  # clamp node 0

geom     = StructureGeometry(nodes, conn, props, constraints)
analysis = FEAAnalysis(geom)

# Load: 2 kN downward at the tip
lc = LoadCase(geom)
lc.add_nodal_load(2, 0, -2000, 0)   # node 2
lc.set_load_factor(1.0)

# Linear analysis
results = analysis.run_linear_analysis(lc)
stress  = analysis.perform_stress_analysis()

print(f"Max disp: {results['max_displacement'] * 1e3:.3f} mm")
print(f"Min FOS:  {stress['FOS']['min_FOS']:.2f}")

# Visualise
import matplotlib.pyplot as plt
from feaopt import StructurePlotter
from feaopt.structure_plotter import draw_beams_stress

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

plotter = StructurePlotter(geom, analysis)

# Deformed shape
deformed = analysis.get_deformed_coordinates(scale=50)
plotter.draw_structure(ax1, color=(0.85, 0.85, 0.85))
plotter.draw_structure(ax1, nodes=deformed, color=(0.2, 0.4, 0.8))
plotter.plot_supports(ax1)
ax1.set_aspect('equal'); ax1.grid(True)
ax1.set_title('Deformed Shape (50x)')

# Stress distribution
vm = stress['stresses']['von_mises']
stress_mpa = np.nanmax(np.abs(vm), axis=1) / 1e6
sm = draw_beams_stress(ax2, geom.nodes, geom.connectivity, geom.thickness, stress_mpa)
plotter.plot_supports(ax2)
plt.colorbar(sm, ax=ax2).set_label('von Mises [MPa]')
ax2.set_aspect('equal'); ax2.grid(True)
ax2.set_title('Stress Distribution')

plt.tight_layout()
plt.show()

This example demonstrates that no file import is necessary — structures can be defined entirely in code by specifying node coordinates and connectivity directly.