oxiphysics_io/

calculix_format.rs

#![allow(clippy::type_complexity)]
// Copyright 2026 COOLJAPAN OU (Team KitaSan)
// SPDX-License-Identifier: Apache-2.0

//! CalculiX `.inp` file format I/O.
//!
//! CalculiX uses an input format that is approximately 95% compatible with
//! the Abaqus `.inp` format.  This module provides a dedicated reader and
//! writer that operate on the generic [`FeMesh`] type from
//! [`crate::finite_element_io`], supporting nodes, elements, materials,
//! boundary conditions, concentrated loads, and analysis steps.

use std::collections::HashMap;
use std::fmt::Write as FmtWrite;
use std::fs;
use std::io::{self, BufRead};

use crate::Error as IoError;
use crate::finite_element_io::{
    AnalysisStep, DirichletBc, FeElement, FeElementType, FeMesh, FeNode, LinearElasticMaterial,
    NodalForce,
};

// ---------------------------------------------------------------------------
// Element type mapping (CalculiX <-> FeElementType)
// ---------------------------------------------------------------------------

/// Map a [`FeElementType`] to its CalculiX string representation.
fn fe_type_to_calculix(et: &FeElementType) -> &'static str {
    match et {
        FeElementType::Tet4 => "C3D4",
        FeElementType::Tet10 => "C3D10",
        FeElementType::Hex8 => "C3D8",
        FeElementType::Hex20 => "C3D20",
        FeElementType::Tri3 => "S3",
        FeElementType::Tri6 => "S6",
        FeElementType::Quad4 => "S4",
        FeElementType::Quad8 => "S8",
        FeElementType::Line2 => "T3D2",
        FeElementType::Unknown(_) => "UNKNOWN",
    }
}

/// Map a CalculiX element-type string to [`FeElementType`] (case-insensitive).
fn calculix_type_to_fe(s: &str) -> FeElementType {
    match s.trim().to_uppercase().as_str() {
        "C3D4" => FeElementType::Tet4,
        "C3D10" => FeElementType::Tet10,
        "C3D8" => FeElementType::Hex8,
        "C3D20" => FeElementType::Hex20,
        "S3" => FeElementType::Tri3,
        "S6" => FeElementType::Tri6,
        "S4" => FeElementType::Quad4,
        "S8" => FeElementType::Quad8,
        "T3D2" => FeElementType::Line2,
        other => FeElementType::Unknown(other.to_string()),
    }
}

// ---------------------------------------------------------------------------
// CalculixWriter
// ---------------------------------------------------------------------------

/// Writes an [`FeMesh`] to a CalculiX `.inp` text file.
#[derive(Debug, Clone, Default)]
pub struct CalculixWriter;

impl CalculixWriter {
    /// Create a new writer.
    pub fn new() -> Self {
        Self
    }

    /// Serialize `mesh` into CalculiX INP format and return the string.
    pub fn write_string(&self, mesh: &FeMesh) -> Result<String, IoError> {
        let mut buf = String::new();

        // --- Heading ---
        writeln!(buf, "*HEADING").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
        writeln!(buf, "Generated by OxiPhysics CalculixWriter")
            .map_err(|e| IoError::General(format!("fmt error: {e}")))?;

        // --- Nodes ---
        writeln!(buf, "*NODE").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
        for node in &mesh.nodes {
            writeln!(
                buf,
                "{}, {:.15e}, {:.15e}, {:.15e}",
                node.id, node.coords[0], node.coords[1], node.coords[2]
            )
            .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
        }

        // --- Elements (grouped by type) ---
        let mut types_seen: Vec<FeElementType> = Vec::new();
        for el in &mesh.elements {
            if !types_seen.contains(&el.element_type) {
                types_seen.push(el.element_type.clone());
            }
        }

        for etype in &types_seen {
            let ccx_type = fe_type_to_calculix(etype);
            writeln!(buf, "*ELEMENT, TYPE={ccx_type}")
                .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            for el in mesh.elements.iter().filter(|e| &e.element_type == etype) {
                let ids: Vec<String> = el.connectivity.iter().map(|n| n.to_string()).collect();
                writeln!(buf, "{}, {}", el.id, ids.join(", "))
                    .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            }
        }

        // --- Element sets (one per element type for section assignment) ---
        if !mesh.element_sets.is_empty() {
            for (set_name, ids) in &mesh.element_sets {
                writeln!(buf, "*ELSET, ELSET={set_name}")
                    .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
                let id_strs: Vec<String> = ids.iter().map(|i| i.to_string()).collect();
                // Write in lines of up to 16 ids
                for chunk in id_strs.chunks(16) {
                    writeln!(buf, "{}", chunk.join(", "))
                        .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
                }
            }
        }

        // --- Materials ---
        for mat in &mesh.materials {
            writeln!(buf, "*MATERIAL, NAME={}", mat.name)
                .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            writeln!(buf, "*DENSITY").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            writeln!(buf, "{:.15e}", mat.density)
                .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            writeln!(buf, "*ELASTIC").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            writeln!(
                buf,
                "{:.15e}, {:.15e}",
                mat.young_modulus, mat.poisson_ratio
            )
            .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
        }

        // --- Solid sections ---
        // Generate a SOLID SECTION for each element set that references a material
        // by convention: element set name == material name
        for mat in &mesh.materials {
            if mesh.element_sets.contains_key(&mat.name) {
                writeln!(
                    buf,
                    "*SOLID SECTION, ELSET={}, MATERIAL={}",
                    mat.name, mat.name
                )
                .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            }
        }

        // --- Steps ---
        for step in &mesh.steps {
            writeln!(buf, "*STEP").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            writeln!(buf, "*STATIC").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
            writeln!(
                buf,
                "{:.6e}, {:.6e}",
                step.initial_increment, step.time_period
            )
            .map_err(|e| IoError::General(format!("fmt error: {e}")))?;

            // Boundary conditions
            if !step.bcs.is_empty() {
                writeln!(buf, "*BOUNDARY")
                    .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
                for bc in &step.bcs {
                    writeln!(
                        buf,
                        "{}, {}, {}, {:.15e}",
                        bc.node_id, bc.dof, bc.dof, bc.value
                    )
                    .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
                }
            }

            // Concentrated loads
            if !step.forces.is_empty() {
                writeln!(buf, "*CLOAD").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
                for f in &step.forces {
                    writeln!(buf, "{}, {}, {:.15e}", f.node_id, f.dof, f.magnitude)
                        .map_err(|e| IoError::General(format!("fmt error: {e}")))?;
                }
            }

            writeln!(buf, "*END STEP").map_err(|e| IoError::General(format!("fmt error: {e}")))?;
        }

        Ok(buf)
    }

    /// Serialize `mesh` and write to the file at `path`.
    pub fn write(&self, mesh: &FeMesh, path: &str) -> Result<(), IoError> {
        let content = self.write_string(mesh)?;
        fs::write(path, content).map_err(IoError::Io)
    }
}

// ---------------------------------------------------------------------------
// CalculixReader
// ---------------------------------------------------------------------------

/// Internal parse-state for the keyword-block state machine.
#[derive(Debug, Clone, PartialEq)]
enum CcxBlock {
    /// Outside any recognised block.
    None,
    /// Inside `*HEADING`.
    Heading,
    /// Inside `*NODE`.
    Node,
    /// Inside `*ELEMENT` with a particular type.
    Element(FeElementType),
    /// Inside `*ELASTIC` for a named material.
    Elastic(String),
    /// Inside `*DENSITY` for a named material.
    Density(String),
    /// Inside `*BOUNDARY`.
    Boundary,
    /// Inside `*CLOAD`.
    Cload,
    /// Inside `*STATIC` (time increment line expected).
    Static,
    /// Inside `*ELSET`.
    Elset(String),
    /// Inside `*NSET`.
    Nset(String),
}

/// Parses CalculiX `.inp` files into an [`FeMesh`].
#[derive(Debug, Clone, Default)]
pub struct CalculixReader;

impl CalculixReader {
    /// Create a new reader.
    pub fn new() -> Self {
        Self
    }

    /// Parse a CalculiX `.inp` file from a file path.
    pub fn parse(&self, path: &str) -> Result<FeMesh, IoError> {
        let file = fs::File::open(path).map_err(IoError::Io)?;
        let reader = io::BufReader::new(file);
        let mut lines = Vec::new();
        for line_res in reader.lines() {
            lines.push(line_res.map_err(IoError::Io)?);
        }
        self.parse_lines(&lines)
    }

    /// Parse a CalculiX `.inp` file from a string.
    pub fn parse_string(&self, source: &str) -> Result<FeMesh, IoError> {
        let lines: Vec<String> = source.lines().map(|l| l.to_string()).collect();
        self.parse_lines(&lines)
    }

    /// Core line-by-line parser.
    fn parse_lines(&self, lines: &[String]) -> Result<FeMesh, IoError> {
        let mut mesh = FeMesh::new();
        let mut block = CcxBlock::None;
        // Track current material name for multi-keyword material definition
        let mut current_material: Option<String> = None;
        // Partial material data: name -> (young, poisson, density)
        let mut mat_data: HashMap<String, (Option<f64>, Option<f64>, Option<f64>)> = HashMap::new();
        // Track whether we are inside a *STEP
        let mut in_step = false;

        for raw_line in lines {
            let line = raw_line.trim();

            // Skip empty lines and comments
            if line.is_empty() || line.starts_with("**") {
                continue;
            }

            // Keyword line?
            if line.starts_with('*') {
                let upper = line.to_uppercase();

                if upper.starts_with("*END STEP") {
                    in_step = false;
                    block = CcxBlock::None;
                    continue;
                }

                if upper.starts_with("*HEADING") {
                    block = CcxBlock::Heading;
                    continue;
                }

                if upper.starts_with("*NODE") && !upper.starts_with("*NSET") {
                    block = CcxBlock::Node;
                    continue;
                }

                if upper.starts_with("*ELEMENT") {
                    let etype_str =
                        extract_param(line, "TYPE").unwrap_or_else(|| "UNKNOWN".to_string());
                    block = CcxBlock::Element(calculix_type_to_fe(&etype_str));
                    continue;
                }

                if upper.starts_with("*NSET") {
                    let name = extract_param(line, "NSET").unwrap_or_else(|| "DEFAULT".to_string());
                    block = CcxBlock::Nset(name);
                    continue;
                }

                if upper.starts_with("*ELSET") {
                    let name =
                        extract_param(line, "ELSET").unwrap_or_else(|| "DEFAULT".to_string());
                    block = CcxBlock::Elset(name);
                    continue;
                }

                if upper.starts_with("*MATERIAL") {
                    let name = extract_param(line, "NAME").unwrap_or_else(|| "UNNAMED".to_string());
                    current_material = Some(name.clone());
                    mat_data.entry(name).or_insert((None, None, None));
                    block = CcxBlock::None;
                    continue;
                }

                if upper.starts_with("*ELASTIC") {
                    if let Some(ref mat_name) = current_material {
                        block = CcxBlock::Elastic(mat_name.clone());
                    } else {
                        block = CcxBlock::None;
                    }
                    continue;
                }

                if upper.starts_with("*DENSITY") {
                    if let Some(ref mat_name) = current_material {
                        block = CcxBlock::Density(mat_name.clone());
                    } else {
                        block = CcxBlock::None;
                    }
                    continue;
                }

                if upper.starts_with("*BOUNDARY") {
                    block = CcxBlock::Boundary;
                    continue;
                }

                if upper.starts_with("*CLOAD") {
                    block = CcxBlock::Cload;
                    continue;
                }

                if upper.starts_with("*STEP") {
                    in_step = true;
                    let name = extract_param(line, "NAME").unwrap_or_else(|| "Step-1".to_string());
                    mesh.steps.push(AnalysisStep::new(name));
                    block = CcxBlock::None;
                    continue;
                }

                if upper.starts_with("*STATIC") {
                    block = CcxBlock::Static;
                    continue;
                }

                if upper.starts_with("*SOLID SECTION") {
                    // We just skip this for now (info is in material + elset)
                    block = CcxBlock::None;
                    continue;
                }

                // Any other keyword: reset block
                block = CcxBlock::None;
                continue;
            }

            // Data line -- dispatch based on current block
            match &block {
                CcxBlock::None | CcxBlock::Heading => {
                    // heading text or unrecognised: ignore
                }
                CcxBlock::Node => {
                    if let Some(node) = parse_node_line(line) {
                        mesh.nodes.push(node);
                    }
                }
                CcxBlock::Element(etype) => {
                    if let Some(elem) = parse_element_line(line, etype.clone()) {
                        mesh.elements.push(elem);
                    }
                }
                CcxBlock::Elastic(mat_name) => {
                    let parts: Vec<&str> = line.split(',').collect();
                    if parts.len() >= 2
                        && let (Some(e), Some(nu)) = (
                            parts[0].trim().parse::<f64>().ok(),
                            parts[1].trim().parse::<f64>().ok(),
                        )
                    {
                        let entry = mat_data
                            .entry(mat_name.clone())
                            .or_insert((None, None, None));
                        entry.0 = Some(e);
                        entry.1 = Some(nu);
                    }
                    block = CcxBlock::None;
                }
                CcxBlock::Density(mat_name) => {
                    if let Some(d) = line
                        .split(',')
                        .next()
                        .and_then(|s| s.trim().parse::<f64>().ok())
                    {
                        let entry = mat_data
                            .entry(mat_name.clone())
                            .or_insert((None, None, None));
                        entry.2 = Some(d);
                    }
                    block = CcxBlock::None;
                }
                CcxBlock::Boundary => {
                    if let Some(bc) = parse_boundary_line(line)
                        && in_step
                        && let Some(step) = mesh.steps.last_mut()
                    {
                        step.bcs.push(bc);
                    }
                    // If outside step, we could store globally; for now skip.
                }
                CcxBlock::Cload => {
                    if let Some(force) = parse_cload_line(line)
                        && in_step
                        && let Some(step) = mesh.steps.last_mut()
                    {
                        step.forces.push(force);
                    }
                }
                CcxBlock::Static => {
                    // Parse time increment line: initial_increment, time_period
                    let parts: Vec<&str> = line.split(',').collect();
                    if parts.len() >= 2
                        && let (Some(inc), Some(period)) = (
                            parts[0].trim().parse::<f64>().ok(),
                            parts[1].trim().parse::<f64>().ok(),
                        )
                        && let Some(step) = mesh.steps.last_mut()
                    {
                        step.initial_increment = inc;
                        step.time_period = period;
                    }
                    block = CcxBlock::None;
                }
                CcxBlock::Elset(name) => {
                    let ids: Vec<usize> = line
                        .split(',')
                        .filter_map(|s| s.trim().parse::<usize>().ok())
                        .collect();
                    mesh.element_sets
                        .entry(name.clone())
                        .or_default()
                        .extend(ids);
                }
                CcxBlock::Nset(name) => {
                    let ids: Vec<usize> = line
                        .split(',')
                        .filter_map(|s| s.trim().parse::<usize>().ok())
                        .collect();
                    mesh.node_sets.entry(name.clone()).or_default().extend(ids);
                }
            }
        }

        // Assemble materials from partial data
        for (mat_id, (name, (young, poisson, density))) in (1usize..).zip(mat_data.iter()) {
            let e = young.unwrap_or(0.0);
            let nu = poisson.unwrap_or(0.0);
            let rho = density.unwrap_or(0.0);
            mesh.materials
                .push(LinearElasticMaterial::new(mat_id, name.clone(), e, nu, rho));
        }

        Ok(mesh)
    }
}

// ---------------------------------------------------------------------------
// Helper parse functions
// ---------------------------------------------------------------------------

/// Extract `KEY=VALUE` from a keyword line (case-insensitive key search).
fn extract_param(line: &str, key: &str) -> Option<String> {
    let upper = line.to_uppercase();
    let key_eq = format!("{}=", key.to_uppercase());
    let pos = upper.find(&key_eq)?;
    let rest = &line[pos + key_eq.len()..];
    let end = rest.find(',').unwrap_or(rest.len());
    let val = rest[..end].trim();
    if val.is_empty() {
        None
    } else {
        Some(val.to_string())
    }
}

/// Parse a node data line: `id, x, y, z`.
fn parse_node_line(line: &str) -> Option<FeNode> {
    let parts: Vec<&str> = line.split(',').collect();
    if parts.len() < 4 {
        return None;
    }
    let id: usize = parts[0].trim().parse().ok()?;
    let x: f64 = parts[1].trim().parse().ok()?;
    let y: f64 = parts[2].trim().parse().ok()?;
    let z: f64 = parts[3].trim().parse().ok()?;
    Some(FeNode::new(id, [x, y, z]))
}

/// Parse an element data line: `id, n1, n2, ...`.
fn parse_element_line(line: &str, etype: FeElementType) -> Option<FeElement> {
    let parts: Vec<&str> = line.split(',').collect();
    if parts.len() < 2 {
        return None;
    }
    let id: usize = parts[0].trim().parse().ok()?;
    let connectivity: Vec<usize> = parts[1..]
        .iter()
        .filter_map(|s| s.trim().parse::<usize>().ok())
        .collect();
    if connectivity.is_empty() {
        return None;
    }
    Some(FeElement::new(id, etype, connectivity))
}

/// Parse a `*BOUNDARY` data line: `node_id, dof_start, dof_end, value`.
///
/// CalculiX boundary lines can be:
/// - `node_id, dof_start, dof_end`          (value = 0)
/// - `node_id, dof_start, dof_end, value`
fn parse_boundary_line(line: &str) -> Option<DirichletBc> {
    let parts: Vec<&str> = line.split(',').collect();
    if parts.len() < 3 {
        return None;
    }
    let node_id: usize = parts[0].trim().parse().ok()?;
    let dof_start: u8 = parts[1].trim().parse().ok()?;
    let dof_end: u8 = parts[2].trim().parse().ok()?;
    let value: f64 = if parts.len() >= 4 {
        parts[3].trim().parse().unwrap_or(0.0)
    } else {
        0.0
    };

    // Expand range: create one BC per DOF in [dof_start, dof_end].
    // For simplicity we return only the first; callers needing the full range
    // should use `parse_boundary_lines`.  However, the typical round-trip
    // writes dof_start == dof_end, so this is fine for most cases.
    // To handle the range properly, we return the first DOF and note that
    // the writer always emits dof_start == dof_end.
    if dof_start == dof_end {
        Some(DirichletBc {
            node_id,
            dof: dof_start,
            value,
        })
    } else {
        // Return just the first DOF of the range
        Some(DirichletBc {
            node_id,
            dof: dof_start,
            value,
        })
    }
}

/// Parse a full `*BOUNDARY` data line into potentially multiple BCs
/// (when `dof_start != dof_end`, one [`DirichletBc`] per DOF in the range).
pub fn parse_boundary_line_expanded(line: &str) -> Vec<DirichletBc> {
    let parts: Vec<&str> = line.split(',').collect();
    if parts.len() < 3 {
        return Vec::new();
    }
    let node_id = match parts[0].trim().parse::<usize>() {
        Ok(v) => v,
        Err(_) => return Vec::new(),
    };
    let dof_start = match parts[1].trim().parse::<u8>() {
        Ok(v) => v,
        Err(_) => return Vec::new(),
    };
    let dof_end = match parts[2].trim().parse::<u8>() {
        Ok(v) => v,
        Err(_) => return Vec::new(),
    };
    let value: f64 = if parts.len() >= 4 {
        parts[3].trim().parse().unwrap_or(0.0)
    } else {
        0.0
    };

    let mut bcs = Vec::new();
    let lo = dof_start.min(dof_end);
    let hi = dof_start.max(dof_end);
    for dof in lo..=hi {
        bcs.push(DirichletBc {
            node_id,
            dof,
            value,
        });
    }
    bcs
}

/// Parse a `*CLOAD` data line: `node_id, dof, magnitude`.
fn parse_cload_line(line: &str) -> Option<NodalForce> {
    let parts: Vec<&str> = line.split(',').collect();
    if parts.len() < 3 {
        return None;
    }
    let node_id: usize = parts[0].trim().parse().ok()?;
    let dof: u8 = parts[1].trim().parse().ok()?;
    let magnitude: f64 = parts[2].trim().parse().ok()?;
    Some(NodalForce::new(node_id, dof, magnitude))
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use std::path::Path;

    // ----- Element type mapping -----

    #[test]
    fn test_fe_type_to_calculix_tet4() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Tet4), "C3D4");
    }

    #[test]
    fn test_fe_type_to_calculix_tet10() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Tet10), "C3D10");
    }

    #[test]
    fn test_fe_type_to_calculix_hex8() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Hex8), "C3D8");
    }

    #[test]
    fn test_fe_type_to_calculix_hex20() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Hex20), "C3D20");
    }

    #[test]
    fn test_fe_type_to_calculix_tri3() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Tri3), "S3");
    }

    #[test]
    fn test_fe_type_to_calculix_tri6() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Tri6), "S6");
    }

    #[test]
    fn test_fe_type_to_calculix_quad4() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Quad4), "S4");
    }

    #[test]
    fn test_fe_type_to_calculix_quad8() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Quad8), "S8");
    }

    #[test]
    fn test_fe_type_to_calculix_line2() {
        assert_eq!(fe_type_to_calculix(&FeElementType::Line2), "T3D2");
    }

    #[test]
    fn test_calculix_type_to_fe_roundtrip() {
        let types = vec![
            FeElementType::Tet4,
            FeElementType::Tet10,
            FeElementType::Hex8,
            FeElementType::Hex20,
            FeElementType::Tri3,
            FeElementType::Tri6,
            FeElementType::Quad4,
            FeElementType::Quad8,
            FeElementType::Line2,
        ];
        for t in types {
            let ccx_str = fe_type_to_calculix(&t);
            let back = calculix_type_to_fe(ccx_str);
            assert_eq!(back, t, "roundtrip failed for {ccx_str}");
        }
    }

    #[test]
    fn test_calculix_type_to_fe_case_insensitive() {
        assert_eq!(calculix_type_to_fe("c3d4"), FeElementType::Tet4);
        assert_eq!(calculix_type_to_fe("C3D10"), FeElementType::Tet10);
        assert_eq!(calculix_type_to_fe("s8"), FeElementType::Quad8);
    }

    #[test]
    fn test_calculix_type_to_fe_unknown() {
        match calculix_type_to_fe("FOOBAR") {
            FeElementType::Unknown(s) => assert_eq!(s, "FOOBAR"),
            other => panic!("expected Unknown, got {:?}", other),
        }
    }

    // ----- Helper functions -----

    #[test]
    fn test_extract_param_type() {
        let line = "*ELEMENT, TYPE=C3D8";
        assert_eq!(extract_param(line, "TYPE"), Some("C3D8".to_string()));
    }

    #[test]
    fn test_extract_param_name() {
        let line = "*MATERIAL, NAME=Steel";
        assert_eq!(extract_param(line, "NAME"), Some("Steel".to_string()));
    }

    #[test]
    fn test_extract_param_missing() {
        let line = "*NODE";
        assert_eq!(extract_param(line, "TYPE"), None);
    }

    #[test]
    fn test_parse_node_line_valid() {
        let n = parse_node_line("1, 0.5, 1.0, 2.0");
        assert!(n.is_some());
        let n = n.expect("should parse");
        assert_eq!(n.id, 1);
        assert!((n.coords[0] - 0.5).abs() < 1e-12);
        assert!((n.coords[1] - 1.0).abs() < 1e-12);
        assert!((n.coords[2] - 2.0).abs() < 1e-12);
    }

    #[test]
    fn test_parse_node_line_invalid() {
        assert!(parse_node_line("1, 0.5").is_none());
    }

    #[test]
    fn test_parse_element_line_valid() {
        let el = parse_element_line("1, 10, 20, 30, 40", FeElementType::Tet4);
        assert!(el.is_some());
        let el = el.expect("should parse");
        assert_eq!(el.id, 1);
        assert_eq!(el.connectivity, vec![10, 20, 30, 40]);
    }

    #[test]
    fn test_parse_element_line_too_short() {
        assert!(parse_element_line("1", FeElementType::Hex8).is_none());
    }

    // ----- Boundary condition parsing -----

    #[test]
    fn test_parse_boundary_line_with_value() {
        let bc = parse_boundary_line("5, 1, 1, 0.01");
        assert!(bc.is_some());
        let bc = bc.expect("should parse");
        assert_eq!(bc.node_id, 5);
        assert_eq!(bc.dof, 1);
        assert!((bc.value - 0.01).abs() < 1e-15);
    }

    #[test]
    fn test_parse_boundary_line_zero_value() {
        let bc = parse_boundary_line("10, 2, 2");
        assert!(bc.is_some());
        let bc = bc.expect("should parse");
        assert_eq!(bc.node_id, 10);
        assert_eq!(bc.dof, 2);
        assert!((bc.value).abs() < 1e-15);
    }

    #[test]
    fn test_parse_boundary_line_invalid() {
        assert!(parse_boundary_line("bad, data").is_none());
    }

    #[test]
    fn test_parse_boundary_line_expanded_range() {
        let bcs = parse_boundary_line_expanded("1, 1, 3, 0.0");
        assert_eq!(bcs.len(), 3);
        assert_eq!(bcs[0].dof, 1);
        assert_eq!(bcs[1].dof, 2);
        assert_eq!(bcs[2].dof, 3);
    }

    // ----- CLOAD parsing -----

    #[test]
    fn test_parse_cload_line_valid() {
        let f = parse_cload_line("7, 2, -1000.0");
        assert!(f.is_some());
        let f = f.expect("should parse");
        assert_eq!(f.node_id, 7);
        assert_eq!(f.dof, 2);
        assert!((f.magnitude - (-1000.0)).abs() < 1e-10);
    }

    #[test]
    fn test_parse_cload_line_invalid() {
        assert!(parse_cload_line("7, 2").is_none());
    }

    // ----- Round-trip: write -> read -> compare -----

    fn sample_mesh() -> FeMesh {
        let mut mesh = FeMesh::new();
        mesh.nodes = vec![
            FeNode::new(1, [0.0, 0.0, 0.0]),
            FeNode::new(2, [1.0, 0.0, 0.0]),
            FeNode::new(3, [0.0, 1.0, 0.0]),
            FeNode::new(4, [0.0, 0.0, 1.0]),
        ];
        mesh.elements = vec![FeElement::new(1, FeElementType::Tet4, vec![1, 2, 3, 4])];
        mesh.materials = vec![LinearElasticMaterial::new(1, "Steel", 210e9, 0.3, 7800.0)];

        let mut step = AnalysisStep::new("LoadStep");
        step.bcs.push(DirichletBc::fixed(1, 1));
        step.bcs.push(DirichletBc::fixed(1, 2));
        step.bcs.push(DirichletBc::fixed(1, 3));
        step.forces.push(NodalForce::new(4, 2, -1000.0));
        step.initial_increment = 0.1;
        step.time_period = 1.0;
        mesh.steps.push(step);

        mesh
    }

    fn tmp_path(name: &str) -> String {
        let dir = std::env::temp_dir();
        dir.join(name).to_string_lossy().to_string()
    }

    #[test]
    fn test_write_creates_file() {
        let path = tmp_path("oxiphysics_ccx_write.inp");
        let mesh = sample_mesh();
        CalculixWriter::new()
            .write(&mesh, &path)
            .expect("write failed");
        assert!(Path::new(&path).exists());
    }

    #[test]
    fn test_roundtrip_node_count() {
        let path = tmp_path("oxiphysics_ccx_rt_nodes.inp");
        let mesh = sample_mesh();
        CalculixWriter::new()
            .write(&mesh, &path)
            .expect("write failed");
        let parsed = CalculixReader::new().parse(&path).expect("parse failed");
        assert_eq!(parsed.nodes.len(), mesh.nodes.len());
    }

    #[test]
    fn test_roundtrip_node_ids() {
        let path = tmp_path("oxiphysics_ccx_rt_nids.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        let ids: Vec<usize> = parsed.nodes.iter().map(|n| n.id).collect();
        assert_eq!(ids, vec![1, 2, 3, 4]);
    }

    #[test]
    fn test_roundtrip_node_coords() {
        let path = tmp_path("oxiphysics_ccx_rt_coords.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert!((parsed.nodes[0].coords[0]).abs() < 1e-10);
        assert!((parsed.nodes[1].coords[0] - 1.0).abs() < 1e-10);
        assert!((parsed.nodes[2].coords[1] - 1.0).abs() < 1e-10);
        assert!((parsed.nodes[3].coords[2] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_roundtrip_element_count() {
        let path = tmp_path("oxiphysics_ccx_rt_elems.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.elements.len(), 1);
    }

    #[test]
    fn test_roundtrip_element_type() {
        let path = tmp_path("oxiphysics_ccx_rt_etype.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.elements[0].element_type, FeElementType::Tet4);
    }

    #[test]
    fn test_roundtrip_element_connectivity() {
        let path = tmp_path("oxiphysics_ccx_rt_econn.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.elements[0].connectivity, vec![1, 2, 3, 4]);
    }

    #[test]
    fn test_roundtrip_material() {
        let path = tmp_path("oxiphysics_ccx_rt_mat.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.materials.len(), 1);
        let mat = &parsed.materials[0];
        assert_eq!(mat.name, "Steel");
        assert!((mat.young_modulus - 210e9).abs() < 1.0);
        assert!((mat.poisson_ratio - 0.3).abs() < 1e-10);
        assert!((mat.density - 7800.0).abs() < 1e-6);
    }

    #[test]
    fn test_roundtrip_boundary_conditions() {
        let path = tmp_path("oxiphysics_ccx_rt_bc.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.steps.len(), 1);
        assert_eq!(parsed.steps[0].bcs.len(), 3);
        assert_eq!(parsed.steps[0].bcs[0].node_id, 1);
        assert_eq!(parsed.steps[0].bcs[0].dof, 1);
        assert!((parsed.steps[0].bcs[0].value).abs() < 1e-15);
    }

    #[test]
    fn test_roundtrip_cload() {
        let path = tmp_path("oxiphysics_ccx_rt_cload.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.steps[0].forces.len(), 1);
        let f = &parsed.steps[0].forces[0];
        assert_eq!(f.node_id, 4);
        assert_eq!(f.dof, 2);
        assert!((f.magnitude - (-1000.0)).abs() < 1e-6);
    }

    #[test]
    fn test_roundtrip_step_timing() {
        let path = tmp_path("oxiphysics_ccx_rt_step.inp");
        let mesh = sample_mesh();
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        let step = &parsed.steps[0];
        assert!((step.initial_increment - 0.1).abs() < 1e-6);
        assert!((step.time_period - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_roundtrip_multiple_element_types() {
        let path = tmp_path("oxiphysics_ccx_rt_multi.inp");
        let mut mesh = FeMesh::new();
        mesh.nodes = vec![
            FeNode::new(1, [0.0, 0.0, 0.0]),
            FeNode::new(2, [1.0, 0.0, 0.0]),
            FeNode::new(3, [0.0, 1.0, 0.0]),
            FeNode::new(4, [0.0, 0.0, 1.0]),
            FeNode::new(5, [1.0, 1.0, 0.0]),
        ];
        mesh.elements = vec![
            FeElement::new(1, FeElementType::Tet4, vec![1, 2, 3, 4]),
            FeElement::new(2, FeElementType::Line2, vec![4, 5]),
        ];
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.elements.len(), 2);
        assert_eq!(parsed.elements[0].element_type, FeElementType::Tet4);
        assert_eq!(parsed.elements[1].element_type, FeElementType::Line2);
    }

    #[test]
    fn test_write_string_contains_heading() {
        let mesh = sample_mesh();
        let content = CalculixWriter::new()
            .write_string(&mesh)
            .expect("write_string");
        assert!(content.contains("*HEADING"));
    }

    #[test]
    fn test_write_string_contains_node() {
        let mesh = sample_mesh();
        let content = CalculixWriter::new()
            .write_string(&mesh)
            .expect("write_string");
        assert!(content.contains("*NODE"));
    }

    #[test]
    fn test_write_string_contains_element() {
        let mesh = sample_mesh();
        let content = CalculixWriter::new()
            .write_string(&mesh)
            .expect("write_string");
        assert!(content.contains("*ELEMENT"));
    }

    #[test]
    fn test_write_string_contains_material() {
        let mesh = sample_mesh();
        let content = CalculixWriter::new()
            .write_string(&mesh)
            .expect("write_string");
        assert!(content.contains("*MATERIAL"));
        assert!(content.contains("*DENSITY"));
        assert!(content.contains("*ELASTIC"));
    }

    #[test]
    fn test_write_string_contains_step() {
        let mesh = sample_mesh();
        let content = CalculixWriter::new()
            .write_string(&mesh)
            .expect("write_string");
        assert!(content.contains("*STEP"));
        assert!(content.contains("*STATIC"));
        assert!(content.contains("*END STEP"));
    }

    #[test]
    fn test_write_string_contains_boundary() {
        let mesh = sample_mesh();
        let content = CalculixWriter::new()
            .write_string(&mesh)
            .expect("write_string");
        assert!(content.contains("*BOUNDARY"));
    }

    #[test]
    fn test_write_string_contains_cload() {
        let mesh = sample_mesh();
        let content = CalculixWriter::new()
            .write_string(&mesh)
            .expect("write_string");
        assert!(content.contains("*CLOAD"));
    }

    #[test]
    fn test_parse_empty_file() {
        let reader = CalculixReader::new();
        let mesh = reader.parse_string("** empty file\n").expect("parse");
        assert!(mesh.nodes.is_empty());
        assert!(mesh.elements.is_empty());
    }

    #[test]
    fn test_parse_missing_file() {
        let result = CalculixReader::new().parse("/tmp/does_not_exist_oxiphysics_ccx.inp");
        assert!(result.is_err());
    }

    #[test]
    fn test_roundtrip_precision() {
        let path = tmp_path("oxiphysics_ccx_rt_prec.inp");
        let mut mesh = FeMesh::new();
        mesh.nodes.push(FeNode::new(
            1,
            [1.23456789012345, -9.87654321098765, 2.89793238462643],
        ));
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        let c = parsed.nodes[0].coords;
        assert!((c[0] - 1.23456789012345).abs() < 1e-10);
        assert!((c[1] - (-9.87654321098765)).abs() < 1e-10);
        assert!((c[2] - 2.89793238462643).abs() < 1e-10);
    }

    #[test]
    fn test_large_mesh_roundtrip() {
        let path = tmp_path("oxiphysics_ccx_rt_large.inp");
        let mut mesh = FeMesh::new();
        for i in 1..=100 {
            mesh.nodes.push(FeNode::new(i, [i as f64, 0.0, 0.0]));
        }
        for i in 0..24usize {
            mesh.elements.push(FeElement::new(
                i + 1,
                FeElementType::Tet4,
                vec![
                    i * 4 % 97 + 1,
                    i * 4 % 97 + 2,
                    i * 4 % 97 + 3,
                    i * 4 % 97 + 4,
                ],
            ));
        }
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.nodes.len(), 100);
        assert_eq!(parsed.elements.len(), 24);
    }

    #[test]
    fn test_hex8_roundtrip() {
        let path = tmp_path("oxiphysics_ccx_rt_hex8.inp");
        let mut mesh = FeMesh::new();
        for i in 1..=8 {
            mesh.nodes.push(FeNode::new(i, [i as f64, 0.0, 0.0]));
        }
        mesh.elements.push(FeElement::new(
            1,
            FeElementType::Hex8,
            vec![1, 2, 3, 4, 5, 6, 7, 8],
        ));
        CalculixWriter::new().write(&mesh, &path).expect("write");
        let parsed = CalculixReader::new().parse(&path).expect("parse");
        assert_eq!(parsed.elements[0].element_type, FeElementType::Hex8);
        assert_eq!(
            parsed.elements[0].connectivity,
            vec![1, 2, 3, 4, 5, 6, 7, 8]
        );
    }

    #[test]
    fn test_parse_string_with_nset_and_elset() {
        let inp = "\
*NODE
1, 0.0, 0.0, 0.0
2, 1.0, 0.0, 0.0
*ELEMENT, TYPE=T3D2
1, 1, 2
*NSET, NSET=FIXED
1
*ELSET, ELSET=ALL
1
";
        let reader = CalculixReader::new();
        let mesh = reader.parse_string(inp).expect("parse");
        assert_eq!(mesh.nodes.len(), 2);
        assert_eq!(mesh.elements.len(), 1);
        assert_eq!(mesh.node_sets.get("FIXED"), Some(&vec![1]));
        assert_eq!(mesh.element_sets.get("ALL"), Some(&vec![1]));
    }

    #[test]
    fn test_parse_boundary_prescribed() {
        let inp = "\
*NODE
1, 0.0, 0.0, 0.0
*STEP
*STATIC
0.01, 1.0
*BOUNDARY
1, 1, 1, 0.05
*END STEP
";
        let reader = CalculixReader::new();
        let mesh = reader.parse_string(inp).expect("parse");
        assert_eq!(mesh.steps.len(), 1);
        assert_eq!(mesh.steps[0].bcs.len(), 1);
        assert_eq!(mesh.steps[0].bcs[0].node_id, 1);
        assert_eq!(mesh.steps[0].bcs[0].dof, 1);
        assert!((mesh.steps[0].bcs[0].value - 0.05).abs() < 1e-15);
    }
}