Layopt

Layopt module.

calc_eq_matrix_b

calc_eq_matrix_b(nodal_coords, c_n, dof)

Calculate equilibrium matrix B.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
c_n	NDArray	Active members.	required
dof	NDArray	Degrees of freedom.	required

Returns:

Type	Description
coo_matrix	Equilibrium matrix B.

Source code in src/layopt/layopt.py

def calc_eq_matrix_b(
    nodal_coords: npt.NDArray, c_n: npt.NDArray, dof: npt.NDArray
) -> sparse.coo_matrix:
    """
    Calculate equilibrium matrix B.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    c_n : npt.NDArray
        Active members.
    dof : npt.NDArray
        Degrees of freedom.

    Returns
    -------
    sparse.coo_matrix
        Equilibrium matrix B.
    """
    try:
        m, n1, n2 = len(c_n), c_n[:, 0].astype(int), c_n[:, 1].astype(int)
    except TypeError as e:
        msg = "Missing 'c_n'"
        raise TypeError(msg) from e

    try:
        length, x, y = (
            c_n[:, 2],
            nodal_coords[n2, 0] - nodal_coords[n1, 0],
            nodal_coords[n2, 1] - nodal_coords[n1, 1],
        )
    except IndexError as e:
        msg = f"{nodal_coords.shape=}, expected (2,{c_n.shape[1]})"
        raise IndexError(msg) from e
    except TypeError as e:
        msg = "Missing 'nodal_coords'"
        raise TypeError(msg) from e

    try:
        d0, d1, d2, d3 = dof[n1 * 2], dof[n1 * 2 + 1], dof[n2 * 2], dof[n2 * 2 + 1]
    except IndexError as e:
        msg = f"{dof.shape=}, expected ({(c_n.shape[0],)})"
        raise IndexError(msg) from e
    except TypeError as e:
        msg = "Missing 'dof'"
        raise TypeError(msg) from e

    s = np.concatenate(
        (-x / length * d0, -y / length * d1, x / length * d2, y / length * d3)
    )
    row_id = np.concatenate((n1 * 2, n1 * 2 + 1, n2 * 2, n2 * 2 + 1))
    col_id = np.concatenate((np.arange(m), np.arange(m), np.arange(m), np.arange(m)))
    return sparse.coo_matrix((s, (row_id, col_id)), shape=(len(nodal_coords) * 2, m))

solve

solve(
    nodal_coords,
    c_n,
    f,
    dof,
    stress_tensile,
    stress_compressive,
    joint_cost,
)

Solve linear programming problem with given connections and pattern load cases.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
c_n	NDArray	Active members.	required
f	list	Load cases.	required
dof	NDArray	Degrees of freedom.	required
stress_tensile	float	Tensile stress limit.	required
stress_compressive	float	Compressive stress limit.	required
joint_cost	float	Joint cost.	required

Returns:

Type	Description
tuple[float, NDArray, list, list]	A tuple consisting of volume (the volume of the solved problem), area (member areas), forces (member forces) and deflections (virtual deflections at degrees of freedom).

Source code in src/layopt/layopt.py

def solve(nodal_coords, c_n, f, dof, stress_tensile, stress_compressive, joint_cost):
    """
    Solve linear programming problem with given connections and pattern load cases.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    c_n : npt.NDArray
        Active members.
    f : list
        Load cases.
    dof : npt.NDArray
        Degrees of freedom.
    stress_tensile : float
        Tensile stress limit.
    stress_compressive : float
        Compressive stress limit.
    joint_cost : float
        Joint cost.

    Returns
    -------
    tuple[float, npt.NDArray, list, list]
        A tuple consisting of ``volume`` (the volume of the solved problem),
        ``area`` (member areas), ``forces`` (member forces) and ``deflections``
        (virtual deflections at degrees of freedom).
    """
    member_cost = [col[2] + joint_cost for col in c_n]
    eq_matrix_b = calc_eq_matrix_b(nodal_coords, c_n, dof)
    q, eqn = [], []
    k = 0  # (damage+load) case number
    with mosek.Model() as model:
        a = model.variable("a", len(c_n), mosek.Domain.greaterThan(0.0))
        eq_matrix_b = mosek.Matrix.sparse(
            eq_matrix_b.shape[0],
            eq_matrix_b.shape[1],
            eq_matrix_b.row,
            eq_matrix_b.col,
            eq_matrix_b.data,
        )
        model.objective(mosek.ObjectiveSense.Minimize, mosek.Expr.dot(member_cost, a))
        for fk in f:
            qi = model.variable(len(c_n))
            q.append(qi)
            eqni = model.constraint(
                mosek.Expr.sub(mosek.Expr.mul(eq_matrix_b, q[k]), fk * dof),
                mosek.Domain.equalsTo(0),
            )
            eqn.append(eqni)
            model.constraint(
                mosek.Expr.sub(mosek.Expr.mul(stress_compressive, a), q[k]),
                mosek.Domain.greaterThan(0),
            )
            model.constraint(
                mosek.Expr.sub(mosek.Expr.mul(-stress_tensile, a), q[k]),
                mosek.Domain.lessThan(0),
            )
            k += 1
        model.setSolverParam("optimizer", "intpnt")
        model.setSolverParam("intpntBasis", "never")
        model.acceptedSolutionStatus(mosek.AccSolutionStatus.Anything)
        # M.setLogHandler(stdout)
        model.solve()
        vol = model.primalObjValue()
        q = [np.array(qi.level()) for qi in q]
        a = a.level()
        u = [-np.array(eqnk.dual()) for eqnk in eqn]
        if vol == 0:
            u = [ui * 10000 for ui in u]

    return vol, a, q, u

stop_violation

stop_violation(
    nodal_coords,
    potential_members,
    dof,
    stress_tensile,
    stress_compressive,
    deflections,
    joint_cost,
)

Check for dual violation and add new members.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
potential_members	NDArray	A list of all possible members.	required
dof	NDArray	Degrees of freedom.	required
stress_tensile	float	Tensile stress limit.	required
stress_compressive	float	Compressive stress limit.	required
deflections	list	Virtual deflections at degrees of freedom.	required
joint_cost	float	Joint cost.	required

Returns:

Type	Description
int	Number of members added.

Source code in src/layopt/layopt.py

def stop_violation(
    nodal_coords: npt.NDArray,
    potential_members: npt.NDArray,
    dof: npt.NDArray,
    stress_tensile: float,
    stress_compressive: float,
    deflections: list,
    joint_cost: float,
):
    """
    Check for dual violation and add new members.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    potential_members : npt.NDArray
        A list of all possible members.
    dof : npt.NDArray
        Degrees of freedom.
    stress_tensile : float
        Tensile stress limit.
    stress_compressive : float
        Compressive stress limit.
    deflections : list
        Virtual deflections at degrees of freedom.
    joint_cost : float
        Joint cost.

    Returns
    -------
    int
        Number of members added.
    """
    lst = np.where(potential_members[:, 3] == False)[0]  # pylint: disable=singleton-comparison
    c_n = potential_members[lst]
    member_cost = c_n[:, 2] + joint_cost
    eq_matrix_b = calc_eq_matrix_b(nodal_coords, c_n, dof).tocsc()
    y = np.zeros(len(c_n))
    for uk in deflections:
        yk = np.multiply(
            eq_matrix_b.transpose().dot(uk) / member_cost,
            np.array([[stress_tensile], [-stress_compressive]]),
        )
        y += np.amax(yk, axis=0)
    vio_c_n = np.where(y > 1.000)[0]
    vio_sort = np.flipud(np.argsort(y[vio_c_n]))
    num = ceil(0.1 * (len(potential_members) - len(c_n)))  # size of existing problem
    for i in range(min(num, len(vio_sort))):
        potential_members[lst[vio_c_n[vio_sort[i]]]][3] = True  # set member as active
    return min(num, len(vio_sort))

plot_truss

plot_truss(
    nodal_coords,
    c_n,
    areas,
    forces,
    threshold,
    title,
    update=True,
    all_cases=False,
)

Visualise truss.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
c_n	NDArray	Active members.	required
areas	NDArray	Member areas.	required
forces	list	Member forces.	required
threshold	float	Minimum allowable member area.	required
title	str	Title for plot, typically includes the filter level expressed as a percentage.	required
update	bool	Enable interactive mode (default=True).	True
all_cases	bool	Plot all load cases individually (default=False).	False

Source code in src/layopt/layopt.py

def plot_truss(
    nodal_coords: npt.NDArray,
    c_n: npt.NDArray,
    areas: list,
    forces: list,
    threshold: float,
    title: str,
    update: bool = True,
    all_cases: bool = False,
):
    """
    Visualise truss.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    c_n : npt.NDArray
        Active members.
    areas : npt.NDArray
        Member areas.
    forces : list
        Member forces.
    threshold : float
        Minimum allowable member area.
    title : str
        Title for plot, typically includes the filter level expressed as a percentage.
    update : bool
        Enable interactive mode (default=``True``).
    all_cases : bool
        Plot all load cases individually (default=``False``).
    """
    if all_cases:
        plot_all_cases(
            nodal_coords=nodal_coords,
            c_n=c_n,
            areas=areas,
            forces=forces,
            threshold=threshold,
            stress_tensile=title,
        )
    else:
        fig = plt.figure()
        ax = fig.subplots()
        if update:
            plt.ion()
        else:
            plt.ioff()
        ax.axis("off")
        ax.axis("equal")
        ax.set_title(title)
        bar_thickness = 0.4  # bar thickness scale
        for i in [i for i in range(len(areas)) if areas[i] >= threshold]:
            if len(forces) > 1:  # multiple LC coloring
                if all(forces[lc][i] >= -0.001 for lc in range(len(forces))):
                    color = "r"
                elif all(forces[lc][i] <= 0.001 for lc in range(len(forces))):
                    color = "b"
                else:
                    color = "tab:gray"
            else:  # single LC coloring (black = no load, dark = less load)
                color = (
                    min(max(forces[0][i] / areas[i], 0), 1),
                    0,
                    min(max(-forces[0][i] / areas[i], 0), 1),
                )
            pos = nodal_coords[c_n[i, [0, 1]].astype(int), :]
            ax.plot(
                pos[:, 0],
                pos[:, 1],
                color=color,
                linewidth=areas[i] * bar_thickness,
                solid_capstyle="round",
            )
        if update:
            plt.pause(0.01)
        else:
            plt.show()

plot_all_cases

plot_all_cases(
    nodal_coords,
    c_n,
    areas,
    forces,
    threshold,
    stress_tensile,
)

Visualise truss, loop through all load cases and plot individually.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
c_n	NDArray	Active members.	required
areas	NDArray	Member areas.	required
forces	list	Member forces.	required
threshold	float	Minimum allowable member area.	required
stress_tensile	str	Stress tensile in string format as a percentage.	required

Source code in src/layopt/layopt.py

def plot_all_cases(
    nodal_coords: npt.NDArray,
    c_n: npt.NDArray,
    areas: npt.NDArray,
    forces: list,
    threshold: float,
    stress_tensile: str,
) -> None:
    """
    Visualise truss, loop through all load cases and plot individually.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    c_n : npt.NDArray
        Active members.
    areas : npt.NDArray
        Member areas.
    forces : list
        Member forces.
    threshold : float
        Minimum allowable member area.
    stress_tensile : str
        Stress tensile in string format as a percentage.
    """
    n_cases = len(forces)
    for k in range(n_cases):
        plot_truss(
            nodal_coords=nodal_coords,
            c_n=c_n,
            areas=areas,
            forces=[forces[k]],
            threshold=threshold,
            title=stress_tensile + " case " + str(k),
            update=True,
            all_cases=False,
        )

make_pattern_loads

make_pattern_loads(
    nodal_coords,
    loaded_points,
    load_large=50.0,
    load_small=5.0,
    load_direction=(0.0, -1.0),
)

Generate all 2^n combinations of large/small loads at each load point.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
loaded_points	NDArray	Load points.	required
load_large	float	Large load to apply at each load point (default=50).	50.0
load_small	float	Small load to apply at each load point (default=5).	5.0
load_direction	tuple[float, float]	Load direction (default=(0,-1)).	(0.0, -1.0)

Returns:

Type	Description
tuple[list, list, list]	A tuple consisting of all_patterns (all load cases), base_load (base load case) and pattern_descriptions (description of each load case using L for large or S for small at each load point).

Source code in src/layopt/layopt.py

def make_pattern_loads(
    nodal_coords: npt.NDArray,
    loaded_points: npt.NDArray,
    load_large: float = 50.0,
    load_small: float = 5.0,
    load_direction: tuple[float, float] = (0.0, -1.0),
):
    """
    Generate all 2^n combinations of large/small loads at each load point.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    loaded_points : npt.NDArray
        Load points.
    load_large : float
        Large load to apply at each load point (default=``50``).
    load_small : float
        Small load to apply at each load point (default=``5``).
    load_direction : tuple[float, float]
        Load direction (default=``(0,-1)``).

    Returns
    -------
    tuple[list, list, list]
        A tuple consisting of ``all_patterns`` (all load cases), ``base_load``
        (base load case) and ``pattern_descriptions`` (description of each load
        case using ``L`` for large or ``S`` for small at each load point).
    """
    if loaded_points is None or not isinstance(loaded_points, np.ndarray):
        msg = f"'loaded_points' is not a numpy array : {type(loaded_points)=}"
        raise TypeError(msg)
    try:
        assert loaded_points.shape[1] >= 1, IndexError(
            f"Need at least one load point : {loaded_points.shape=}"
        )
    except IndexError as e:
        msg = f"Need at least one load point : {loaded_points.shape}"
        raise IndexError(msg) from e

    # Find node indices for each load point
    load_node_indices = []
    for loaded_point in loaded_points:
        dists = np.linalg.norm(nodal_coords - np.array(loaded_point), axis=1)
        load_node_indices.append(np.argmin(dists))

    # ns-rse 2026-03-17 : Could maybe use a dictionary here to link patterns and descriptions?
    all_patterns = []
    pattern_descriptions = []

    # Generate all 2^n combinations
    # First combo (all loadLarge) is base case
    for combo in itertools.product([load_large, load_small], repeat=len(loaded_points)):
        fk = np.zeros(len(nodal_coords) * 2)
        desc = []
        for pt_idx, magnitude in enumerate(combo):
            node = load_node_indices[pt_idx]
            fk[node * 2] += magnitude * load_direction[0]
            fk[node * 2 + 1] += magnitude * load_direction[1]
            desc.append(f"pt{pt_idx}={'L' if magnitude == load_large else 'S'}")

        all_patterns.append(fk)
        pattern_descriptions.append(", ".join(desc))

    # ns-rse 2026-03-17 : Return directly as part of tuple
    base_load = all_patterns[0]  # First pattern = all large loads
    logger.info(
        f"Total patterns for {len(loaded_points)} load point(s) : {len(all_patterns)}"
    )
    logger.info(f"Base case (all large) : {pattern_descriptions[0]}")

    return all_patterns, base_load, pattern_descriptions

stop_primal_violation_residual

stop_primal_violation_residual(
    nodal_coords,
    c_n,
    forces,
    all_patterns,
    active_load_cases,
    dof,
)

Check for primal violation (equilibrium constraint violation) and add new load cases.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
c_n	NDArray	Active members.	required
forces	list	Member forces.	required
all_patterns	list	All load cases.	required
active_load_cases	list	For each load case, bool set to True if active, False otherwise.	required
dof	NDArray	Degrees of freedom.	required

Returns:

Type	Description
bool	True if converged and no load cases added.

Source code in src/layopt/layopt.py

def stop_primal_violation_residual(
    nodal_coords: npt.NDArray,
    c_n: npt.NDArray,
    forces: npt.NDArray,
    all_patterns: list,
    active_load_cases: list,
    dof: npt.NDArray,
) -> bool:
    """
    Check for primal violation (equilibrium constraint violation) and add new load cases.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    c_n : npt.NDArray
        Active members.
    forces : list
        Member forces.
    all_patterns : list
        All load cases.
    active_load_cases : list
        For each load case, bool set to ``True`` if active, ``False`` otherwise.
    dof : npt.NDArray
        Degrees of freedom.

    Returns
    -------
    bool
        True if converged and no load cases added.
    """
    tol = 1e-5
    eq_matrix_b = calc_eq_matrix_b(nodal_coords, c_n, dof).tocsc()

    total_violation = np.zeros(len(all_patterns))

    # loop through all (active and inactive) pattern load cases
    for k, _ in enumerate(all_patterns):
        if active_load_cases[k] == 1:
            continue  # skip active cases

        fk_dof = all_patterns[k] * dof
        residuals = [
            np.linalg.norm(eq_matrix_b.dot(force) - fk_dof) for force in forces
        ]

        # find min of residuals over active pattern load cases
        total_violation[k] = min(residuals)

    violated = (
        total_violation > tol
    )  # true if there are violated cases need to be added
    n_to_add = max(1, ceil(len(all_patterns) / 10))  # limit on num to add

    if any(violated):
        # ns-rse 2026-03-17 : extract sorting violations to its own function
        # Sort by violation severity
        by_violation = sorted(
            [i for i in range(len(total_violation)) if total_violation[i] > tol],
            key=lambda k: total_violation[k],
            reverse=True,
        )

        if len(by_violation) == 0:
            return True

        # Active most violated load pattern
        active_load_cases[by_violation[0]] = 1
        violations_added_this_iter = [by_violation[0]]
        by_violation.pop(0)

        # Add distinct cases
        # distinct if violated load pattern vector (after normalisation)
        # is not parallel to added load pattern vector (after normalisation)
        # (with current loading everything should be distinct?)
        for _ in range(n_to_add - 1):
            if len(by_violation) == 0:
                break
            added_case = False
            for k in by_violation:
                fk = all_patterns[k]
                fk_norm = fk / (np.linalg.norm(fk) + 1e-12)
                distinct = True
                for j in violations_added_this_iter:
                    fj = all_patterns[j]
                    fj_norm = fj / (np.linalg.norm(fj) + 1e-12)
                    if np.dot(fk_norm, fj_norm) > 0.99:
                        distinct = False
                        break
                if distinct:
                    active_load_cases[k] = 1
                    violations_added_this_iter.append(k)
                    by_violation.remove(k)
                    added_case = True
                    break
            if not added_case:
                break

        return False  # cases added, keep going
    return True  # converged, terminate

stop_primal_violation_pattern

stop_primal_violation_pattern(
    nodal_coords,
    c_n,
    areas,
    all_patterns,
    active_load_cases,
    dof,
    stress_tensile,
    stress_compressive,
)

Check for primal violation (load factor structural analysis) and add new load cases.

Parameters:

Name	Type	Description	Default
nodal_coords	NDArray	Nodal coordinates.	required
c_n	NDArray	Active members.	required
areas	list[float]	Member areas.	required
all_patterns	list	All load cases.	required
active_load_cases	list	For each load case, bool set to True if active, False otherwise.	required
dof	NDArray	Degrees of freedom.	required
stress_tensile	float	Tensile stress limit.	required
stress_compressive	float	Compressive stress limit.	required

Returns:

Type	Description
bool	True if converged and no load cases added.

Source code in src/layopt/layopt.py

def stop_primal_violation_pattern(
    nodal_coords: npt.NDArray,
    c_n: npt.NDArray,
    areas: list[float],
    all_patterns: list,
    active_load_cases: list,
    dof: npt.NDArray,
    stress_tensile: float,
    stress_compressive: float,
) -> bool:
    """
    Check for primal violation (load factor structural analysis) and add new load cases.

    Parameters
    ----------
    nodal_coords : npt.NDArray
        Nodal coordinates.
    c_n : npt.NDArray
        Active members.
    areas : list[float]
        Member areas.
    all_patterns : list
        All load cases.
    active_load_cases : list
        For each load case, bool set to True if active, False otherwise.
    dof : npt.NDArray
        Degrees of freedom.
    stress_tensile : float
        Tensile stress limit.
    stress_compressive : float
        Compressive stress limit.

    Returns
    -------
    bool
        ``True`` if converged and no load cases added.
    """
    tol = 0.99  # lambda must be >= 1 to be considered feasible

    eq_matrix_b = calc_eq_matrix_b(nodal_coords, c_n, dof)
    load_factors = np.ones(len(all_patterns))  # lambda=1 for active cases

    # loop through all (active and inactive) pattern load cases
    for k, _ in enumerate(all_patterns):
        if active_load_cases[k] == 1:
            continue  # skip active cases

        fk_dof = all_patterns[k] * dof

        # Solve LP: maximize lambda subject to B*q = lambda*f, -sigma*a <= q <= sigma*a
        with mosek.Model() as model:
            # Variables
            q_var = model.variable("q", len(c_n))
            lambda_var = model.variable("lambda", 1, mosek.Domain.greaterThan(0.0))

            # Objective: maximize lambda
            model.objective(mosek.ObjectiveSense.Maximize, lambda_var)

            # Constraint 1: B*q = lambda*f
            b_mosek = mosek.Matrix.sparse(
                eq_matrix_b.shape[0],
                eq_matrix_b.shape[1],
                eq_matrix_b.row,
                eq_matrix_b.col,
                eq_matrix_b.data,
            )
            model.constraint(
                mosek.Expr.sub(
                    mosek.Expr.mul(b_mosek, q_var),
                    mosek.Expr.mul(fk_dof, lambda_var.index(0)),
                ),
                mosek.Domain.equalsTo(0),
            )

            # Constraint 2: q <= sigma_t * a (tension limit)
            model.constraint(
                mosek.Expr.sub(q_var, stress_tensile * areas), mosek.Domain.lessThan(0)
            )

            # Constraint 3: q >= -sigma_c * a (compression limit)
            model.constraint(
                mosek.Expr.sub(q_var, -stress_compressive * areas),
                mosek.Domain.greaterThan(0),
            )

            # Solve
            model.setSolverParam("optimizer", "intpnt")
            model.acceptedSolutionStatus(mosek.AccSolutionStatus.Anything)
            model.solve()
            load_factors[k] = lambda_var.level()[0]

    # Violation: load factor < 1 (with tolerance)
    violated = load_factors < tol
    n_to_add = max(1, ceil(len(all_patterns) / 10))

    if any(violated):  # pylint: disable=too-many-nested-blocks
        # ns-rse 2026-03-17 : extract sorting violations to its own function
        # Sort by severity: lowest load factor = most violated
        by_violation = sorted(
            [i for i in range(len(load_factors)) if violated[i]],
            key=lambda k: load_factors[k],
        )

        if len(by_violation) == 0:
            return True

        # Add most violated (lowest lambda)
        most_violated_id = by_violation[0]
        active_load_cases[most_violated_id] = 1
        logger.info(
            f"  Adding most violated pattern {by_violation[0]}: lambda={load_factors[most_violated_id]:.3f}"
        )
        violations_added_this_iter = [by_violation[0]]
        by_violation.pop(0)

        # Add distinct cases
        # distinct if load factor is +/-10% of added load factor
        # and if violated load pattern vector (after normalisation)
        # is not parallel to added load pattern vector (after normalisation)
        # (with current loading no load pattern vectors should be parallel?)
        for _ in range(n_to_add - 1):
            if len(by_violation) == 0:
                break
            added_case = False
            for k in by_violation:
                # check if load case k has a significantly different load factor
                # from all load cases added this iteration
                distinct = True
                for j in violations_added_this_iter:
                    # check if both load factors are approx 0, only add one if so
                    if load_factors[k] < 0.01 and load_factors[j] < 0.01:
                        distinct = False
                        break
                    # if added load factor is 0 but other violated ones aren't,
                    # other violated cases are distinct
                    if load_factors[j] < 0.01:
                        continue
                    # check ratio of load factors if neither approx 0
                    lambda_ratio = load_factors[k] / (load_factors[j] + 1e-12)
                    if 0.9 < lambda_ratio < 1.1:  # Within 10% of each other
                        distinct = False
                        break
                if distinct:
                    fk = all_patterns[k]
                    fk_norm = fk / (np.linalg.norm(fk) + 1e-12)
                    distinct = True
                    for j in violations_added_this_iter:
                        fj = all_patterns[j]
                        fj_norm = fj / (np.linalg.norm(fj) + 1e-12)
                        if np.dot(fk_norm, fj_norm) > 0.99:
                            distinct = False
                            break
                if distinct:
                    active_load_cases[k] = 1
                    logger.info(
                        f"  Adding {len(violations_added_this_iter) + 1} distinct pattern {k}: lambda={load_factors[k]:.3f}"
                    )
                    violations_added_this_iter.append(k)
                    by_violation.remove(k)
                    added_case = True
                    break
            if not added_case:
                break

        return False  # cases added, keep going
    return True  # converged, terminate

trussopt

trussopt(
    filter_level=None,
    width=1.0,
    height=1.0,
    stress_tensile=1.0,
    stress_compressive=1.0,
    joint_cost=0.0,
    loaded_points=None,
    load_direction=(0.0, -1.0),
    load_large=50.0,
    load_small=5.0,
    max_length=1000.0,
    support_points=None,
    primal_method="load_factor",
    problem_name="None",
    notes="",
)

Main function, perform adaptive member adding procedure with multiple load cases.

Parameters:

Name	Type	Description	Default
filter_level	float	Levels to filter on.	None
width	float	Width of structure.	1.0
height	float	Height of structure.	1.0
stress_tensile	float	Tensile stress limit.	1.0
stress_compressive	float	Compressive stress limit.	1.0
joint_cost	float	Joint cost.	0.0
loaded_points	list	Load points (default=[]).	None
load_direction	list	Load direction (default=(0,-1)).	(0.0, -1.0)
load_large	float	Large load to apply at each load point (default=50).	50.0
load_small	float	Small load to apply at each load point (default=5).	5.0
max_length	float	Maximum member length.	1000.0
support_points	NDArray	Support points (default=[]).	None
primal_method	str	Primal violation method (default='load_factor').	'load_factor'
problem_name	str	Name of problem to solve (default=None).	'None'
notes	str	Notes (default='').	''

Returns:

Type	Description
tuple[float, NDArray, DataFrame, float]	A tuple consisting of volume (the final volume of the solved problem) and area (final member areas of the solved problem), a dataframe of results and the filter_level.

Source code in src/layopt/layopt.py

def trussopt(
    filter_level: float | None = None,
    width: float = 1.0,
    height: float = 1.0,
    stress_tensile: float = 1.0,
    stress_compressive: float = 1.0,
    joint_cost: float = 0.0,
    loaded_points: list | None = None,
    # ns-rse 2026-03-17 : val implies single value but its a list, perhaps load_range? dict perhaps
    load_direction: tuple[float, float] = (0.0, -1.0),
    load_large: float = 50.0,
    load_small: float = 5.0,
    max_length: float = 1000.0,
    # ns-rse 2026-03-17 : Set type hint and default to None
    support_points: npt.NDArray | None = None,
    primal_method: str = "load_factor",
    problem_name: str = "None",
    # save_to_csv: bool = True,
    # csv_filename: str = "pattern_loading_results.csv",
    notes: str = "",
) -> tuple[float, npt.NDArray, pd.DataFrame, float]:
    """
    Main function, perform adaptive member adding procedure with multiple load cases.

    Parameters
    ----------
    filter_level : float
        Levels to filter on.
    width : float
        Width of structure.
    height : float
        Height of structure.
    stress_tensile : float
        Tensile stress limit.
    stress_compressive : float
        Compressive stress limit.
    joint_cost : float
        Joint cost.
    loaded_points : list
        Load points (default=[]).
    load_direction : list
        Load direction (default=``(0,-1)``).
    load_large : float
        Large load to apply at each load point (default=50).
    load_small : float
        Small load to apply at each load point (default=5).
    max_length : float
        Maximum member length.
    support_points : npt.NDArray
        Support points (default=[]).
    primal_method : str
        Primal violation method (default='load_factor').
    problem_name : str
        Name of problem to solve (default=``None``).
    notes : str
        Notes (default='').

    Returns
    -------
    tuple[float, npt.NDArray, pd.DataFrame, float]
        A tuple consisting of ``volume`` (the final volume of the solved problem)
        and ``area`` (final member areas of the solved problem), a dataframe of results and the ``filter_level``.
    """
    setup_start = time.process_time()

    # Make domain
    poly = Polygon([(0, 0), (width, 0), (width, height), (0, height)])
    convex = poly.convex_hull.area == poly.area
    logger.debug(f"Domain created, convex? : {convex=}")

    # Make nodes
    xv, yv = np.meshgrid(range(width + 1), range(height + 1))
    points = [Point(xv.flat[i], yv.flat[i]) for i in range(xv.size)]
    logger.debug(f"Points created : {len(points)=}")
    nodal_coords = np.array([[pt.x, pt.y] for pt in points if poly.intersects(pt)])
    logger.debug(f"Node coordinates :\n{nodal_coords=}")
    dof = np.ones((len(nodal_coords), 2))
    logger.debug(f"Degrees of Freedom : {dof=}")

    # Default load point
    if loaded_points is None:
        loaded_points = [[width, height // 2]]
        logger.info("Loaded points not provided, calculated as : {loaded_points=}")
    # support conditions
    for i, node in enumerate(nodal_coords):
        if support_points.size == 0:
            if node[0] == 0:
                dof[i, :] = [0, 0]  # Support nodes with x=0
        else:
            dof[i, :] = (
                [0, 0]
                if any((node == point).all() for point in support_points)
                else [1, 1]
            )
    dof = np.array(dof).flatten()

    # Generate all pattern loads
    # ns-rse 2026-03-17 : Unused arguments but may combine all_patterns and pattern_descriptions to dict
    # all_patterns, base_load, pattern_descriptions = make_pattern_loads(
    all_patterns, _, _ = make_pattern_loads(
        nodal_coords, loaded_points, load_large, load_small, load_direction
    )

    # Create the 'ground structure'
    potential_members = []
    for i, j in itertools.combinations(range(len(nodal_coords)), 2):
        dx, dy = (
            abs(nodal_coords[i][0] - nodal_coords[j][0]),
            abs(nodal_coords[i][1] - nodal_coords[j][1]),
        )
        length = np.sqrt(dx**2 + dy**2)
        # Remove overlapping members, or members longer than maxLength
        if (length < max_length and gcd(int(dx), int(dy)) == 1) or joint_cost != 0:
            seg = [] if convex else LineString([nodal_coords[i], nodal_coords[j]])
            if convex or poly.contains(seg) or poly.boundary.contains(seg):
                potential_members.append([i, j, length, False])
    potential_members = np.array(potential_members)

    # Create the active members
    # DualAdaptivity = True
    # start_len = 1.5 if DualAdaptivity else 10000
    # for pm in [p for p in PML if p[2] <= start_len]:  # Activate short members (adaptive)
    # ns-rse 2026-03-16 : DualAdaptivity is always 'True'
    for pm in [
        p for p in potential_members if p[2] <= 1.5
    ]:  # Activate short members (adaptive)
        pm[3] = True

    #### Primal adaptivity: start with base load case only ####
    primal_adaptivity = True
    # if PrimalAdaptivity:
    #     activeLoadCases = np.zeros(len(allPatterns), dtype=int)
    #     activeLoadCases[0] = 1  # Base case = all large loads
    # # does below make sense here?
    # else:
    #     activeLoadCases = np.ones(len(allPatterns), dtype=int)
    if primal_method in {"residual", "load_factor"}:
        primal_adaptivity = True
        active_load_cases = np.zeros(len(all_patterns), dtype=int)
        active_load_cases[0] = 1  # Base case = all large loads
    else:
        primal_adaptivity = False
        active_load_cases = np.ones(len(all_patterns), dtype=int)

    setup_end = time.process_time()
    logger.info(f"Setup took {setup_end - setup_start!s}")
    logger.info(f"    Nodes               : {len(nodal_coords)}")
    logger.info(f"    Members             : {len(potential_members)}")
    logger.info(f"    Total load patterns : {len(all_patterns)}")

    vol = 1e9  # arbitrary large number to initialise
    # Start the 'member adding' loop
    for itr in range(1, 100):
        last_volume = vol
        # Get active members/parts of matrices
        c_n = potential_members[potential_members[:, 3] == True]  # pylint: disable=singleton-comparison

        # Get active pattern loads
        f_active = [
            all_patterns[k]
            for k in range(len(all_patterns))
            if active_load_cases[k] == 1
        ]

        # solve current reduced problem
        vol, a, q, u = solve(
            nodal_coords,
            c_n,
            f_active,
            dof,
            stress_tensile,
            stress_compressive,
            joint_cost,
        )
        # We need to solve once so that we have valid values for `a` which we then filter based on `fitler_level[s]`
        # (rename to `filter_level` but need to check first if that is what we want to parallelise on or if it is
        # `primal_method`).
        # if filter_level != 1.0:
        #     # max_a = max(a)
        #     # filter_val = filter_level * max_a
        #     # keep = [a_value > filter_val for a_value in a]
        #     keep = [_a > (filter_level * max(a)) for _a in a]
        #     kept = c_n[keep]
        # c_n = [_a for _a in a if _a > filter_level * max(a)]

        # output
        if isinf(vol):
            logger.error("Infeasible problem detected")
            return [], [], []
        n_active = int(np.sum(active_load_cases))
        # ns-rse 2026-03-23 : Could this perhaps be debugging?
        logger.info(
            f"Iteration: {itr}, vol: {vol}, mems: {len(c_n)} active load cases:{n_active}/{len(all_patterns)}"
        )
        # plot interim solutions (slow)
        # plotTruss(nodal_coords, c_n, a, q, max(a) * 1e-2, "Itr:" + str(itr), extraPlot = activeDamageDef)

        # inner loop - adding of members based on dual violation
        # still need PMLcache? currently unused
        # PMLcache = np.copy(PML[:,3])
        n_added = stop_violation(
            nodal_coords,
            potential_members,
            dof,
            stress_tensile,
            stress_compressive,
            u,
            joint_cost,
        )
        if not (0.99 * last_volume) < vol < (1.0001 * last_volume):
            continue  # small vol decrease = member adding close to convergence

        # outer loop - adding of pattern load cases based on primal violation
        # if stopPrimalViolationPattern(nodal_coords, c_n, a, all_patterns, active_load_cases, dof, st, sc):
        #     if numAdded > 0: # only fully terminate when no members violate
        #         continue
        #     else:
        #         break

        if primal_adaptivity:
            if primal_method == "residual":
                # Use equilibrium residual check
                converged = stop_primal_violation_residual(
                    nodal_coords, c_n, q, all_patterns, active_load_cases, dof
                )
            elif primal_method == "load_factor":
                # Use load factor LP
                converged = stop_primal_violation_pattern(
                    nodal_coords,
                    c_n,
                    a,
                    all_patterns,
                    active_load_cases,
                    dof,
                    stress_tensile,
                    stress_compressive,
                )
            # ns-rse 2026-03-17 : leaves scope for 'converged' to not be assigned if `primal_method` never matches

            if not converged:  # pylint: disable=possibly-used-before-assignment
                continue  # Cases added, keep iterating
            if n_added > 0:
                continue  # No cases added but members added
            break  # Both converged
        # No primal adaptivity - just check member convergence
        if n_added == 0:
            break  # Converged

    final_vol = vol
    logger.info(f"Volume: {final_vol}")
    solve_end = time.process_time()
    logger.info("Solve took " + str(solve_end - setup_end))
    logger.info(
        f"Active patterns: {int(np.sum(active_load_cases))}/{len(all_patterns)}"
    )
    # Build dictionary of results
    results = {
        "timestamp": get_date_time(),
        "problem_name": problem_name or f"w{width}_h{height}_n{len(loaded_points)}",
        "filter_level": filter_level,
        "width": width,
        "height": height,
        "n_load_points": len(loaded_points),
        "n_patterns_total": len(all_patterns),
        "n_patterns_active": int(np.sum(active_load_cases)),
        "load_large": load_large,
        "load_small": load_small,
        "iterations": itr,
        "final_volume": None,
        "n_members_final": len(c_n),
        "n_nodes": len(nodal_coords),
        "n_ground_structure": len(potential_members),
        "cpu_time_setup": setup_end - setup_start,
        "cpu_time_solve": solve_end - setup_end,
        # 'cpu_time_total': total_cpu_time,
        # 'wall_time_total': total_wall_time,
        "primal_method": primal_method,
        "notes": notes,
    }

    # Plot results
    # plotTruss(nodal_coords, c_n, a, q, max(a)*1e-2, "Final", update=False, allCases=True)

    ## Filter
    # if filter_levels:
    #     results["final_vol"] = {}
    #     for multiplier in filter_levels:
    #         max_a = max(a)
    #         filter_val = multiplier * max_a
    #         keep = [a_value > filter_val for a_value in a]
    #         kept = c_n[keep]
    #         vol, filer_a, filter_q, u = solve(
    #             nodal_coords,
    #             kept,
    #             f_active,
    #             dof,
    #             stress_tensile,
    #             stress_compressive,
    #             joint_cost,
    #         )
    #         if vol > 0:
    #             logger.info(
    #                 f"filtered volume {vol} with filter at {100 * multiplier}% gives {len(filer_a)!s} members"
    #             )
    #             plot_truss(
    #                 nodal_coords=nodal_coords,
    #                 c_n=kept,
    #                 areas=filer_a,
    #                 forces=filter_q,
    #                 threshold=max(a) * 1e-3,
    #                 title="Filtered " + str(100 * multiplier) + "%",
    #                 update=False,
    #                 all_cases=False,
    #             )
    #         results["final_vol"][multiplier] = vol

    #     logger.info(f"Plotting took {time.process_time() - solve_end!s}")
    #     save_results_to_csv(results, csv_filename)
    #     return vol, a, results, None

    results["final_volume"] = final_vol
    # Save results to CSV
    # if save_to_csv:
    #     # ns-rse 2026-03-16 - inefficient to write to CSV, build dictionary/dataframe in memory and write to disk on
    #     #                     completion (as we may end up paralllelising processing)
    #     save_results_to_csv(results, csv_filename)
    return vol, a, dict_to_df(results), filter_level

save_results_to_csv

save_results_to_csv(
    results, filename="pattern_loading_results.csv"
)

Save optimization results to CSV file.

Creates file with header if it doesn't exist, otherwise appends.

Parameters:

Name	Type	Description	Default
results	dict	Dictionary containing results to save. Should include keys: - timestamp - problem_name - width, height - n_load_points - n_patterns_total - n_patterns_active - load_large - load_small - iterations - final_volume - n_members_final - n_nodes - n_ground_structure - cpu_time_setup - cpu_time_solve - primal_method - notes	required
filename	str	CSV filename (default=pattern_loading_results.csv).	'pattern_loading_results.csv'

Source code in src/layopt/layopt.py

def save_results_to_csv(
    results: dict[str, str | int | float], filename="pattern_loading_results.csv"
):
    """
    Save optimization results to CSV file.

    Creates file with header if it doesn't exist, otherwise appends.

    Parameters
    ----------
    results : dict
        Dictionary containing results to save. Should include keys:
        - timestamp
        - problem_name
        - width, height
        - n_load_points
        - n_patterns_total
        - n_patterns_active
        - load_large
        - load_small
        - iterations
        - final_volume
        - n_members_final
        - n_nodes
        - n_ground_structure
        - cpu_time_setup
        - cpu_time_solve
        - primal_method
        - notes
    filename : str
        CSV filename (default=``pattern_loading_results.csv``).
    """
    file_exists = Path(filename).is_file()

    # Define column order
    fieldnames = [
        "timestamp",
        "problem_name",
        "width",
        "height",
        "n_load_points",
        "n_patterns_total",
        "n_patterns_active",
        "load_large",
        "load_small",
        "iterations",
        "final_volume",
        "n_members_final",
        "n_nodes",
        "n_ground_structure",
        "cpu_time_setup",
        "cpu_time_solve",
        # 'cpu_time_total',
        # 'wall_time_total',
        "primal_method",
        "notes",
    ]

    with Path(filename).open("a", newline="", encoding="utf-8") as csvfile:
        writer = csv.DictWriter(csvfile, fieldnames=fieldnames)

        # Write header if new file
        if not file_exists:
            writer.writeheader()

        # Write data
        writer.writerow(results)

    logger.info(f"\nResults saved to {filename}")