Structural Analysis (CalculiX)

frond.fem_solver

CalculiX FEM Simulation Builder.

Reads a Gmsh-generated mesh .inp file and appends the complete CalculiX simulation deck: material definitions, shell sections, orientations, MPC equations, boundary conditions, loads, and analysis step.

build_ccx_deck(mesh_inp: str, web_properties: Dict[int, dict], segments: list = None, skin_thickness: float = 0.003, web_thickness_override: float = None, tip_load: float = -1000.0, point_loads: List[Dict[str, float]] = None, material: dict = None, output_path: str = None, solver: str = None, binary_output: bool = False, mapped_aero_forces: Dict[int, np.ndarray] = None, skin_elset_name: str = 'WingSkin') -> str

Build a complete CalculiX simulation deck from a Gmsh mesh .inp file.

Parameters:

Name	Type	Description	Default
mesh_inp	str	Path to the Gmsh-generated .inp mesh file (already has NODE, ELEMENT, *ELSET).	required
web_properties	dict	Per-web property dict from build_brep_webs(). Keys are web indices (0, 1, ...), values are dicts with 'thickness', 'level', 'id'.	required
segments	list of Seg	The fractal segments, used to compute per-web fiber orientation vectors.	None
skin_thickness	float	Shell thickness for the wing skin (default: 3mm).	0.003
web_thickness_override	float	If set, overrides all web thicknesses with this value.	None
tip_load	float	Fallback concentrated force if point_loads is not provided.	-1000.0
point_loads	list of dict	List of dictionaries defining point loads at the tip. Format: [{'x_frac': 0.5, 'load': -1000.0}], where x_frac is the chordwise position (0=LE, 1=TE).	None
material	dict	Material properties dict. Defaults to CFRP_T300.	None
output_path	str	Path for the output simulation .inp file. Defaults to replacing the mesh file's extension with '_sim.inp'.	None
solver	str	Equation solver keyword (e.g. 'PARDISO', 'PASTIX'). If None, CalculiX uses its default solver (SPOOLES). Only effective when the ccx binary has been compiled with the corresponding library.	None
binary_output	bool	If True, use *NODE OUTPUT / *ELEMENT OUTPUT for compact mixed binary/ASCII .frd files (~40%% smaller, faster I/O). If False (default), use *NODE FILE / *EL FILE for fully ASCII output compatible with ccx2paraview and other post-processing tools. Set to True only when VTU conversion is not needed.	False

Returns:

Type	Description
str	Path to the generated simulation .inp file.

Source code in src/frond/fem_solver.py

def build_ccx_deck(
    mesh_inp: str,
    web_properties: Dict[int, dict],
    segments: list = None,
    skin_thickness: float = 0.003,
    web_thickness_override: float = None,
    tip_load: float = -1000.0,
    point_loads: List[Dict[str, float]] = None,
    material: dict = None,
    output_path: str = None,
    solver: str = None,
    binary_output: bool = False,
    mapped_aero_forces: Dict[int, np.ndarray] = None,
    skin_elset_name: str = 'WingSkin',
) -> str:
    """
    Build a complete CalculiX simulation deck from a Gmsh mesh .inp file.

    Parameters
    ----------
    mesh_inp : str
        Path to the Gmsh-generated .inp mesh file (already has *NODE, *ELEMENT, *ELSET).
    web_properties : dict
        Per-web property dict from build_brep_webs(). Keys are web indices (0, 1, ...),
        values are dicts with 'thickness', 'level', 'id'.
    segments : list of Seg, optional
        The fractal segments, used to compute per-web fiber orientation vectors.
    skin_thickness : float
        Shell thickness for the wing skin (default: 3mm).
    web_thickness_override : float, optional
        If set, overrides all web thicknesses with this value.
    tip_load : float
        Fallback concentrated force if point_loads is not provided.
    point_loads : list of dict, optional
        List of dictionaries defining point loads at the tip.
        Format: [{'x_frac': 0.5, 'load': -1000.0}], where x_frac is the chordwise position (0=LE, 1=TE).
    material : dict, optional
        Material properties dict. Defaults to CFRP_T300.
    output_path : str, optional
        Path for the output simulation .inp file. Defaults to replacing the
        mesh file's extension with '_sim.inp'.
    solver : str, optional
        Equation solver keyword (e.g. ``'PARDISO'``, ``'PASTIX'``).  If None,
        CalculiX uses its default solver (SPOOLES).  Only effective when the
        ``ccx`` binary has been compiled with the corresponding library.
    binary_output : bool
        If True, use ``*NODE OUTPUT`` / ``*ELEMENT OUTPUT`` for compact
        mixed binary/ASCII ``.frd`` files (~40%% smaller, faster I/O).
        If False (default), use ``*NODE FILE`` / ``*EL FILE`` for fully
        ASCII output compatible with ``ccx2paraview`` and other
        post-processing tools.  Set to True only when VTU conversion
        is not needed.

    Returns
    -------
    str
        Path to the generated simulation .inp file.
    """
    if material is None:
        material = CFRP_T300

    if output_path is None:
        output_path = mesh_inp.replace('.inp', '_sim.inp')

    # ── Parse the mesh ──
    nodes = _parse_inp_nodes(mesh_inp)
    elements = _parse_inp_elements(mesh_inp)
    elsets = _parse_inp_elsets(mesh_inp)

    # ── Identify coordinate-based node sets ──
    ys = [c[1] for c in nodes.values()]
    y_min, y_max = min(ys), max(ys)
    tol_y = 1e-3  # 1mm tolerance for coordinate matching

    root_nodes = sorted(
        [nid for nid, (x, y, z) in nodes.items() if abs(y - y_min) < tol_y])
    tip_nodes = sorted(
        [nid for nid, (x, y, z) in nodes.items() if abs(y - y_max) < tol_y])

    # Process point loads at the tip
    active_loads = []
    if point_loads is not None:
        for p in point_loads:
            frac = p.get('x_frac', 0.5)
            val = p.get('load', tip_load)
            active_loads.append({'x_frac': frac, 'load': val})
    else:
        active_loads.append({'x_frac': 0.5, 'load': tip_load})

    applied_loads = []
    if tip_nodes:
        tip_xs = [nodes[n][0] for n in tip_nodes]
        min_x, max_x = min(tip_xs), max(tip_xs)
        dx = max_x - min_x if max_x > min_x else 1.0

        for p in active_loads:
            target_x = min_x + p['x_frac'] * dx
            node_id = min(tip_nodes, key=lambda n: abs(nodes[n][0] - target_x))
            applied_loads.append((node_id, p['load']))
    else:
        # Fallback if no tip nodes found
        fallback_node = max(nodes.keys())
        for p in active_loads:
            applied_loads.append((fallback_node, p['load']))

    # ── Detect Web and Beam Elsets and Boundary Nodes ──
    web_elset_names = sorted([k for k in set(list(elsets.keys()) + list(elements.keys())) if k.startswith('Web_')])
    beam_elset_names = sorted([k for k in set(list(elsets.keys()) + list(elements.keys())) if k.startswith('BEAMS_T_')])

    # Flatten elements dict for quick lookup by ID
    flat_elements = {}
    for elem_list in elements.values():
        for eid, enodes in elem_list:
            flat_elements[eid] = enodes

    # Topologically find all web edges
    web_edges = {}
    for wname in web_elset_names:
        for eid in elsets.get(wname, []):
            if eid in flat_elements:
                enodes = flat_elements[eid]
                n_n = len(enodes)
                for i in range(n_n):
                    edge = tuple(sorted((enodes[i], enodes[(i + 1) % n_n])))
                    web_edges[edge] = web_edges.get(edge, 0) + 1

    # Nodes on edges shared by only 1 web element are on the outer boundary
    boundary_nodes = set()
    for edge, count in web_edges.items():
        if count == 1:
            boundary_nodes.add(edge[0])
            boundary_nodes.add(edge[1])

    # For 2D beams, EVERY beam node must be tied to the skin
    for bname in beam_elset_names:
        eids = set(elsets.get(bname, []))
        if bname in elements:
            eids.update([eid for eid, _ in elements[bname]])
        for eid in eids:
            if eid in flat_elements:
                boundary_nodes.update(flat_elements[eid])

    # Exclude root nodes to prevent SPC/MPC conflicts
    boundary_nodes.difference_update(root_nodes)

    # ── Match Web Boundary Nodes to Skin ──
    skin_nodes = _get_element_nodes_for_elset(skin_elset_name, elements, elsets)
    skin_nodes_list = list(skin_nodes)
    skin_coords = np.array([nodes[n] for n in skin_nodes_list if n in nodes])

    mpc_equations = []
    if len(skin_coords) > 0 and boundary_nodes:
        from scipy.spatial import cKDTree
        skin_tree = cKDTree(skin_coords)

        w_nodes_list = list(boundary_nodes)
        w_coords = np.array([nodes[n] for n in w_nodes_list if n in nodes])

        if len(w_coords) == len(w_nodes_list):
            _, indices = skin_tree.query(w_coords, k=1)
            for i in range(len(w_nodes_list)):
                web_n = w_nodes_list[i]
                skin_n = skin_nodes_list[indices[i]]
                mpc_equations.append((web_n, skin_n))

    # ── Compute per-web orientation vectors ──
    web_orientations = {}
    if segments is not None:
        for web_idx, props in web_properties.items():
            seg_id = props.get('id', web_idx)
            if seg_id < len(segments):
                seg = segments[seg_id]
                direction = seg.p1 - seg.p0
                length = np.linalg.norm(direction)
                if length > 1e-12:
                    direction = direction / length
                    web_orientations[web_idx] = (
                        float(
                            direction[0]), float(
                            direction[1]), 0.0)

    # ── Build the simulation deck ──
    lines = []

    # Read the original mesh content
    with open(mesh_inp, 'r') as f:
        mesh_content = f.read()
    lines.append(mesh_content.rstrip())
    lines.append('')

    # ── Node Sets ──
    lines.append('**')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('** NODE SETS')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('**')

    # Root nodes
    lines.append('*NSET, NSET=NSET_ROOT')
    for k in range(0, len(root_nodes), 10):
        chunk = root_nodes[k:k + 10]
        lines.append(', '.join(str(n) for n in chunk) + ',')

    # Tip nodes
    lines.append('*NSET, NSET=NSET_TIP')
    for k in range(0, len(tip_nodes), 10):
        chunk = tip_nodes[k:k + 10]
        lines.append(', '.join(str(n) for n in chunk) + ',')

    # Tip load nodes
    lines.append('*NSET, NSET=NSET_TIP_LOAD')
    for n_id, _ in applied_loads:
        lines.append(f'{n_id},')

    # ── Material Definition ──
    lines.append('**')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('** MATERIAL')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('**')
    lines.append(f'*MATERIAL, NAME={material["name"]}')
    lines.append('*ELASTIC, TYPE=ENGINEERING CONSTANTS')
    lines.append(f'{material["E1"]}, {material["E2"]}, {material["E3"]}, '
                 f'{material["nu12"]}, {material["nu13"]}, {material["nu23"]}, '
                 f'{material["G12"]}, {material["G13"]},')
    lines.append(f'{material["G23"]}')
    lines.append('*DENSITY')
    lines.append(f'{material["density"]}')

    # ── Orientations ──
    lines.append('**')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('** ORIENTATIONS')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('**')

    lines.append('*ORIENTATION, NAME=OR_SKIN, SYSTEM=RECTANGULAR')
    lines.append('0., 1., 0.,  0., 0., 1.')

    for web_idx in range(len(web_elset_names)):
        if web_idx in web_orientations:
            dx, dy, dz = web_orientations[web_idx]
        else:
            dx, dy, dz = 0., 1., 0.
        lines.append(
            f'*ORIENTATION, NAME=OR_Web_{web_idx}, SYSTEM=RECTANGULAR')
        lines.append(f'{dx}, {dy}, {dz},  0., 0., 1.')

    # ── Shell Sections ──
    lines.append('**')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('** SHELL SECTIONS')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('**')

    lines.append(f'*SHELL SECTION, ELSET={skin_elset_name}, MATERIAL={material["name"]}, '
                 f'ORIENTATION=OR_SKIN')
    lines.append(f'{skin_thickness}')

    for web_idx in range(len(web_elset_names)):
        if web_thickness_override is not None:
            thick = web_thickness_override
        elif web_idx in web_properties:
            thick = web_properties[web_idx].get('thickness', 0.005)
        else:
            thick = 0.005
        lines.append(f'*SHELL SECTION, ELSET=Web_{web_idx}, MATERIAL={material["name"]}, '
                     f'ORIENTATION=OR_Web_{web_idx}')
        lines.append(f'{thick}')

    # ── Beam Sections ──
    if beam_elset_names:
        lines.append('**')
        lines.append('** ═══════════════════════════════════════════')
        lines.append('** BEAM SECTIONS')
        lines.append('** ═══════════════════════════════════════════')
        lines.append('**')
        for bname in beam_elset_names:
            try:
                # Name is like BEAMS_T_0_003000
                thick_str = bname.replace('BEAMS_T_', '')
                # Find where the decimal point should be based on length, or just replace last underscore
                # Actually, in meshing.py we did f"BEAMS_T_{thick:.6f}".replace('.', '_')
                # So if it's 0_003000, we can replace the first underscore with a dot
                parts = thick_str.split('_')
                if len(parts) >= 2:
                    thick = float(parts[0] + '.' + parts[1])
                else:
                    thick = 0.005
            except (ValueError, IndexError):
                thick = 0.005

            lines.append(f'*BEAM SECTION, ELSET={bname}, MATERIAL={material["name"]}, SECTION=CIRC')
            lines.append(f'{thick / 2.0}')  # CIRC takes radius
            lines.append('0., 0., -1.')  # Direction cosines for beam orientation (z-axis downwards)

    # ── MPC Equations (Web to Skin Contact) ──
    if mpc_equations:
        lines.append('**')
        lines.append('** ═══════════════════════════════════════════')
        lines.append('** KINEMATIC COUPLING (MPC EQUATIONS)')
        lines.append('** ═══════════════════════════════════════════')
        lines.append('** Tying web boundary nodes to nearest skin nodes')
        for web_n, skin_n in mpc_equations:
            for dof in range(1, 7):
                lines.append('*EQUATION')
                lines.append('2')
                lines.append(f'{web_n}, {dof}, 1.0, {skin_n}, {dof}, -1.0')

    # ── Boundary Conditions ──
    lines.append('**')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('** BOUNDARY CONDITIONS AND LOADS')
    lines.append('** ═══════════════════════════════════════════')
    lines.append('**')

    lines.append('*STEP')
    if solver:
        lines.append(f'*STATIC, SOLVER={solver.upper()}')
    else:
        lines.append('*STATIC')
    lines.append('1., 1.')
    lines.append('**')

    lines.append('*BOUNDARY')
    lines.append('NSET_ROOT, 1, 6')
    lines.append('**')

    if mapped_aero_forces is not None:
        lines.append('**')
        lines.append('** MAPPED AERODYNAMIC LOADS')
        lines.append('**')
        for n_id, force in mapped_aero_forces.items():
            lines.append('*CLOAD')
            if abs(force[0]) > 1e-6:
                lines.append(f'{n_id}, 1, {force[0]:.8e}'.upper())
            if abs(force[1]) > 1e-6:
                lines.append(f'{n_id}, 2, {force[1]:.8e}'.upper())
            if abs(force[2]) > 1e-6:
                lines.append(f'{n_id}, 3, {force[2]:.8e}'.upper())
        lines.append('**')
    else:
        for n_id, load_val in applied_loads:
            lines.append('*CLOAD')
            lines.append(f'{n_id}, 3, {load_val}')
            lines.append('**')

    if binary_output:
        lines.append('*NODE OUTPUT')
        lines.append('U')
        lines.append('*ELEMENT OUTPUT')
        lines.append('S')
    else:
        lines.append('*NODE FILE')
        lines.append('U')
        lines.append('*EL FILE')
        lines.append('S')
    lines.append('**')
    lines.append('*END STEP')

    with open(output_path, 'w') as f:
        f.write('\n'.join(lines) + '\n')

    print(f'  -> CalculiX simulation deck written: {output_path}')
    print(f'     Root nodes (clamped): {len(root_nodes)}')
    print(f'     Tip nodes: {len(tip_nodes)}')
    for n_id, load_val in applied_loads:
        print(f'     Tip load node: {n_id} ({load_val} N in Z)')
    print(f'     Web-to-Skin MPCs: {len(mpc_equations)}')
    print(f'     Web orientations defined: {len(web_orientations)}')

    return output_path

parse_mesh_for_mapping(mesh_inp: str, elset_name: str = 'WingSkin') -> Tuple[Dict[int, Tuple[float, float, float]], Set[int]]

Parse mesh file to get all nodes and the set of skin node IDs.

Source code in src/frond/fem_solver.py

def parse_mesh_for_mapping(
        mesh_inp: str, elset_name: str = 'WingSkin') -> Tuple[Dict[int, Tuple[float, float, float]], Set[int]]:
    """Parse mesh file to get all nodes and the set of skin node IDs."""
    nodes = _parse_inp_nodes(mesh_inp)
    elements = _parse_inp_elements(mesh_inp)
    elsets = _parse_inp_elsets(mesh_inp)
    skin_nodes = _get_element_nodes_for_elset(elset_name, elements, elsets)
    return nodes, skin_nodes