caexfer 0.1.2

//! Read supported MSH datasets.

use crate::core::{
    Cell, CellKind, CellSet, Dataset, Error, Field, FieldLocation, Mesh, NodeSet, Point, Result,
};
use std::collections::{BTreeMap, BTreeSet};
use std::fmt::Write as _;

/// Mesh projection and metadata that had no native group mapping.
#[derive(Debug, Clone, PartialEq)]
pub struct Projection {
    /// Complete supported geometry and numeric fields.
    pub dataset: Dataset,

    /// Elements with nonzero geometrical or other unmapped 2.2 tags.
    pub tagged_elements: usize,

    /// Physical groups assigned deterministic names because none were stored.
    pub generated_group_names: usize,
}

/// One MSH 2.2 element before regrouping into internal 4.1 blocks.
struct Element22 {
    /// Original element tag.
    id: u64,

    /// Supported linear Gmsh type code.
    code: u32,

    /// Original node tags in Gmsh ordering.
    nodes: Vec<u64>,

    /// First metadata tag, when it is a positive physical group number.
    physical: Option<u64>,
}

/// Attach the MSH-specific diagnostic code to a parsing or writing failure.
pub(super) fn err(message: impl Into<String>) -> Error {
    Error::new("E_MSH", message)
}

/// Parse a required token, preserving its field name in the diagnostic.
fn number<T: std::str::FromStr>(value: Option<&str>, what: &str) -> Result<T> {
    value
        .ok_or_else(|| err(format!("missing {what}")))?
        .parse()
        .map_err(|_| err(format!("invalid {what}")))
}

/// Decode a Gmsh linear element code; reject unsupported element orders.
fn kind(code: u32) -> Result<CellKind> {
    match code {
        1 => Ok(CellKind::Line2),
        2 => Ok(CellKind::Triangle3),
        3 => Ok(CellKind::Quad4),
        4 => Ok(CellKind::Tet4),
        5 => Ok(CellKind::Hex8),
        6 => Ok(CellKind::Wedge6),
        7 => Ok(CellKind::Pyramid5),
        _ => Err(err(format!("unsupported element type {code}"))),
    }
}

/// Return the Gmsh element code and topological dimension for a cell.
pub(super) fn code(kind: CellKind) -> (u32, u32) {
    match kind {
        CellKind::Line2 => (1, 1),
        CellKind::Triangle3 => (2, 2),
        CellKind::Quad4 => (3, 2),
        CellKind::Tet4 => (4, 3),
        CellKind::Hex8 => (5, 3),
        CellKind::Wedge6 => (6, 3),
        CellKind::Pyramid5 => (7, 3),
    }
}

/// Borrow the first named section body, checking that its end marker exists.
fn section<'a>(source: &'a str, name: &str) -> Result<Option<&'a str>> {
    // A missing section is optional to callers; a present truncated one is not.
    let begin = format!("${name}\n");
    let end = format!("$End{name}");
    let Some(start) = source.find(&begin) else {
        return Ok(None);
    };
    let body_start = start + begin.len();
    let stop = source[body_start..]
        .find(&end)
        .ok_or_else(|| err(format!("unterminated ${name}")))?
        + body_start;
    Ok(Some(&source[body_start..stop]))
}

/// Borrow all repeated field sections in source order, rejecting truncation.
fn data_sections<'a>(source: &'a str, name: &str) -> Result<Vec<&'a str>> {
    // NodeData and ElementData may occur once per field and step.
    let mut out = Vec::new();
    let mut rest = source;
    let begin = format!("${name}\n");
    let end = format!("$End{name}");
    while let Some(start) = rest.find(&begin) {
        rest = &rest[start + begin.len()..];
        let stop = rest
            .find(&end)
            .ok_or_else(|| err(format!("unterminated ${name}")))?;
        out.push(&rest[..stop]);
        rest = &rest[stop + end.len()..];
    }
    Ok(out)
}

/// Decode explicit physical group names, retaining spaces inside quotes.
fn physical_names(source: &str) -> Result<BTreeMap<(u8, u64), String>> {
    let Some(body) = section(source, "PhysicalNames")? else {
        return Ok(BTreeMap::new());
    };
    let mut lines = body.lines().map(str::trim).filter(|line| !line.is_empty());
    let count: usize = number(lines.next(), "physical name count")?;
    if count > body.len() / 5 {
        return Err(err("physical name count exceeds input"));
    }
    let mut names = BTreeMap::new();
    for _ in 0..count {
        let line = lines.next().ok_or_else(|| err("missing physical name"))?;
        let mut words = line.split_whitespace();
        let dimension: u8 = number(words.next(), "physical dimension")?;
        let tag: u64 = number(words.next(), "physical tag")?;
        let start = line
            .find('"')
            .ok_or_else(|| err("physical name must be quoted"))?;
        let end = line
            .rfind('"')
            .ok_or_else(|| err("physical name must be quoted"))?;
        if start == end || dimension > 3 || tag == 0 {
            return Err(err("invalid physical name record"));
        }
        let name = line[start + 1..end].to_owned();
        if name.is_empty() || names.insert((dimension, tag), name).is_some() {
            return Err(err("empty or duplicate physical group name"));
        }
    }
    if lines.next().is_some() {
        return Err(err("extra physical name records"));
    }
    Ok(names)
}

/// Read entity physical tags from the structured MSH 4.1 entity section.
fn entity_physical_tags(source: &str) -> Result<BTreeMap<(u8, u64), Vec<u64>>> {
    let Some(body) = section(source, "Entities")? else {
        return Ok(BTreeMap::new());
    };
    let mut lines = body.lines().map(str::trim).filter(|line| !line.is_empty());
    let counts = lines
        .next()
        .ok_or_else(|| err("missing entity counts"))?
        .split_whitespace()
        .map(|word| {
            word.parse::<usize>()
                .map_err(|_| err("invalid entity count"))
        })
        .collect::<Result<Vec<_>>>()?;
    let counts: [usize; 4] = counts
        .try_into()
        .map_err(|_| err("four entity counts required"))?;
    let declared = counts
        .iter()
        .try_fold(0_usize, |total, count| total.checked_add(*count))
        .ok_or_else(|| err("entity count overflows"))?;
    if declared > body.len() / 5 {
        return Err(err("entity count exceeds input"));
    }
    let mut entities = BTreeMap::new();
    for (dimension, count) in counts.into_iter().enumerate() {
        for _ in 0..count {
            let words = lines
                .next()
                .ok_or_else(|| err("missing entity record"))?
                .split_whitespace()
                .collect::<Vec<_>>();
            let prefix = if dimension == 0 { 4 } else { 7 };
            if words.len() <= prefix {
                return Err(err("short entity record"));
            }
            let tag: u64 = number(words.first().copied(), "entity tag")?;
            let physical_count: usize =
                number(words.get(prefix).copied(), "entity physical count")?;
            if physical_count > words.len().saturating_sub(prefix + 1) {
                return Err(err("entity physical count exceeds record"));
            }
            let physical = words[prefix + 1..prefix + 1 + physical_count]
                .iter()
                .map(|word| {
                    word.parse::<u64>()
                        .map_err(|_| err("invalid entity physical tag"))
                })
                .collect::<Result<Vec<_>>>()?;
            let dimension = u8::try_from(dimension).map_err(|_| err("invalid entity dimension"))?;
            if entities.insert((dimension, tag), physical).is_some() {
                return Err(err("duplicate entity tag"));
            }
        }
    }
    if lines.next().is_some() {
        return Err(err("extra entity records"));
    }
    Ok(entities)
}

/// Attach overlapping physical groups to the cells they select.
fn attach_cell_sets(
    mesh: &mut Mesh,
    names: &BTreeMap<(u8, u64), String>,
    groups: BTreeMap<(u8, u64), Vec<u64>>,
) -> usize {
    let mut generated = 0;
    for ((dimension, tag), cell_ids) in groups {
        let name = names.get(&(dimension, tag)).cloned().unwrap_or_else(|| {
            generated += 1;
            format!("physical_{dimension}_{tag}")
        });
        mesh.cell_sets.push(CellSet {
            name,
            dimension,
            cell_ids,
        });
    }
    generated
}

/// Attach physical point groups from dimension-zero node entities.
fn attach_node_sets(
    mesh: &mut Mesh,
    names: &BTreeMap<(u8, u64), String>,
    groups: BTreeMap<u64, Vec<u64>>,
) -> usize {
    let mut generated = 0;
    for (tag, point_ids) in groups {
        let name = names.get(&(0, tag)).cloned().unwrap_or_else(|| {
            generated += 1;
            format!("physical_0_{tag}")
        });
        mesh.node_sets.push(NodeSet { name, point_ids });
    }
    generated
}

/// Read an ASCII MSH 4.1 file with linear elements and numeric data blocks.
///
/// Binary files, parametric nodes, and unknown element types are rejected.
/// Node and element tags remain the original IDs; connectivity becomes indices
/// into [`Mesh::points`].
///
/// # Errors
///
/// Returns an error for malformed sections, declared counts exceeding input,
/// unsupported element types, or incomplete numeric fields.
///
#[allow(clippy::too_many_lines)] // MSH section counts and records are validated in one pass.
fn read_41(source: &str) -> Result<(Dataset, usize)> {
    let names = physical_names(source)?;
    let entity_physical = entity_physical_tags(source)?;
    let mut selected_cells: BTreeMap<(u8, u64), Vec<u64>> = BTreeMap::new();
    let mut selected_nodes: BTreeMap<u64, Vec<u64>> = BTreeMap::new();

    // Gmsh stores node tags before coordinate tuples within each block.
    // Build a tag-to-index map for later element connectivity.
    let mut mesh = Mesh::default();
    let mut node_ids = BTreeMap::new();
    let nodes = section(source, "Nodes")?.ok_or_else(|| err("missing Nodes"))?;
    let mut t = nodes.split_whitespace();
    let blocks: usize = number(t.next(), "node block count")?;
    let total: usize = number(t.next(), "node count")?;
    let _: u64 = number(t.next(), "minimum node tag")?;
    let _: u64 = number(t.next(), "maximum node tag")?;
    if blocks > t.clone().count() / 4 || total > t.clone().count() / 4 {
        return Err(err("node counts exceed available input"));
    }
    for _ in 0..blocks {
        let entity_dim: u8 = number(t.next(), "entity dimension")?;
        let entity_tag: u64 = number(t.next(), "entity tag")?;
        let parametric: u8 = number(t.next(), "parametric flag")?;
        if parametric != 0 {
            return Err(err("parametric nodes are unsupported"));
        }
        let count: usize = number(t.next(), "block node count")?;
        if count > total.saturating_sub(mesh.points.len()) || count > t.clone().count() / 4 {
            return Err(err("block node count exceeds available input"));
        }
        let mut tags = Vec::with_capacity(count);
        for _ in 0..count {
            tags.push(number::<u64>(t.next(), "node tag")?);
        }
        for id in tags {
            let position = [
                number(t.next(), "x")?,
                number(t.next(), "y")?,
                number(t.next(), "z")?,
            ];
            if node_ids.insert(id, mesh.points.len()).is_some() {
                return Err(err("duplicate node tag"));
            }
            mesh.points.push(Point { id, position });
            if entity_dim == 0 {
                if let Some(tags) = entity_physical.get(&(0, entity_tag)) {
                    for &tag in tags {
                        selected_nodes.entry(tag).or_default().push(id);
                    }
                }
            }
        }
    }
    if mesh.points.len() != total {
        return Err(err("node count mismatch"));
    }

    // Element blocks are homogeneous in type and entity dimension.
    if let Some(elements) = section(source, "Elements")? {
        let mut t = elements.split_whitespace();
        let blocks: usize = number(t.next(), "element block count")?;
        let total: usize = number(t.next(), "element count")?;
        let _: u64 = number(t.next(), "minimum element tag")?;
        let _: u64 = number(t.next(), "maximum element tag")?;
        if blocks > t.clone().count() / 4 || total > t.clone().count() / 3 {
            return Err(err("element counts exceed available input"));
        }
        for _ in 0..blocks {
            let dim: u32 = number(t.next(), "entity dimension")?;
            let entity_tag: u64 = number(t.next(), "entity tag")?;
            let type_code: u32 = number(t.next(), "element type")?;
            let cell_kind = kind(type_code)?;
            if dim != code(cell_kind).1 {
                return Err(err("element dimension mismatch"));
            }
            let count: usize = number(t.next(), "block element count")?;
            let tokens_per_element = cell_kind.node_count() + 1;
            if count > total.saturating_sub(mesh.cells.len())
                || count > t.clone().count() / tokens_per_element
            {
                return Err(err("block element count exceeds available input"));
            }
            for _ in 0..count {
                let id = number(t.next(), "element tag")?;
                let mut connectivity = Vec::with_capacity(cell_kind.node_count());
                for _ in 0..cell_kind.node_count() {
                    let node: u64 = number(t.next(), "element node")?;

                    // The mesh model stores point indices, not Gmsh node tags.
                    connectivity.push(
                        *node_ids
                            .get(&node)
                            .ok_or_else(|| err(format!("unknown node {node}")))?,
                    );
                }
                mesh.cells.push(Cell {
                    id,
                    kind: cell_kind,
                    connectivity,
                    property_id: None,
                });
                if let Some(tags) = entity_physical.get(&(
                    u8::try_from(dim).map_err(|_| err("invalid entity dimension"))?,
                    entity_tag,
                )) {
                    for &tag in tags {
                        selected_cells
                            .entry((
                                u8::try_from(dim).map_err(|_| err("invalid entity dimension"))?,
                                tag,
                            ))
                            .or_default()
                            .push(id);
                    }
                }
            }
        }
        if mesh.cells.len() != total {
            return Err(err("element count mismatch"));
        }
    }
    let generated_group_names = attach_cell_sets(&mut mesh, &names, selected_cells)
        + attach_node_sets(&mut mesh, &names, selected_nodes);
    let mut dataset = Dataset {
        mesh,
        fields: Vec::new(),
    };

    // Numeric data blocks identify entities by original tag and can appear
    // in an order different from the geometry arrays.
    for (name, location) in [
        ("NodeData", FieldLocation::Point),
        ("ElementData", FieldLocation::Cell),
    ] {
        for body in data_sections(source, name)? {
            let mut lines = body.lines().map(str::trim).filter(|line| !line.is_empty());
            let strings: usize = number(lines.next(), "string tag count")?;
            if strings != 1 {
                return Err(err("data block requires exactly one name tag"));
            }
            let raw_name = lines.next().ok_or_else(|| err("missing field name"))?;
            if !raw_name.starts_with('"') || !raw_name.ends_with('"') {
                return Err(err("field name must be quoted"));
            }
            let field_name = raw_name[1..raw_name.len() - 1].to_string();
            let reals: usize = number(lines.next(), "real tag count")?;
            let time: Option<f64> = match reals {
                0 => None,
                1 => Some(number(lines.next(), "time")?),
                _ => return Err(err("data block supports at most one time tag")),
            };
            let ints: usize = number(lines.next(), "integer tag count")?;
            if ints != 3 {
                return Err(err("data block requires three integer tags"));
            }
            let step: i64 = number(lines.next(), "step")?;
            let components: usize = number(lines.next(), "component count")?;
            let entries: usize = number(lines.next(), "data entry count")?;
            let remainder = lines.collect::<Vec<_>>().join(" ");
            let mut t = remainder.split_whitespace();
            if components == 0 {
                return Err(err("zero-component field"));
            }
            let values_len = entries
                .checked_mul(components)
                .ok_or_else(|| err("data value count overflows usize"))?;
            if components > source.len() || entries > t.clone().count() / (components + 1) {
                return Err(err("data counts exceed available input"));
            }
            let ids: Vec<u64> = match location {
                FieldLocation::Point => dataset.mesh.points.iter().map(|p| p.id).collect(),
                FieldLocation::Cell => dataset.mesh.cells.iter().map(|c| c.id).collect(),
            };
            if entries != ids.len() {
                return Err(err(
                    "partial data blocks cannot be projected without missing-value semantics",
                ));
            }

            // Reorder each complete block into the mesh's entity-major order.
            let index: BTreeMap<u64, usize> =
                ids.iter().enumerate().map(|(i, id)| (*id, i)).collect();
            let mut values = vec![0.; values_len];
            let mut seen = BTreeSet::new();
            for _ in 0..entries {
                let id: u64 = number(t.next(), "data entity tag")?;
                let at = *index
                    .get(&id)
                    .ok_or_else(|| err("data references an unknown entity"))?;
                if !seen.insert(id) {
                    return Err(err("duplicate data entity"));
                }
                for component in 0..components {
                    values[at * components + component] = number(t.next(), "data value")?;
                }
            }
            dataset.fields.push(Field {
                name: field_name,
                location,
                components: (0..components).map(|i| format!("C{}", i + 1)).collect(),
                values,
                step: Some(step),
                time,
            });
        }
    }
    dataset.validate()?;
    Ok((dataset, generated_group_names))
}

/// Parse MSH 2.2 flat sections and reuse the common data-block decoder.
#[allow(clippy::too_many_lines)] // Flat node/element records each require bounded count checks.
fn read_22(source: &str) -> Result<Projection> {
    let nodes = section(source, "Nodes")?.ok_or_else(|| err("missing Nodes"))?;
    let mut tokens = nodes.split_whitespace();
    let count: usize = number(tokens.next(), "node count")?;
    if count > tokens.clone().count() / 4 {
        return Err(err("MSH 2.2 node count exceeds available input"));
    }
    let mut points = Vec::with_capacity(count);
    for _ in 0..count {
        let id: u64 = number(tokens.next(), "node tag")?;
        let coords: [f64; 3] = [
            number(tokens.next(), "x")?,
            number(tokens.next(), "y")?,
            number(tokens.next(), "z")?,
        ];
        points.push((id, coords));
    }
    if tokens.next().is_some() {
        return Err(err("extra MSH 2.2 node tokens"));
    }

    let elements = section(source, "Elements")?.ok_or_else(|| err("missing Elements"))?;
    let mut tokens = elements.split_whitespace();
    let count: usize = number(tokens.next(), "element count")?;
    if count > tokens.clone().count() / 4 {
        return Err(err("MSH 2.2 element count exceeds available input"));
    }
    let mut parsed = Vec::with_capacity(count);
    let mut tagged_elements = 0usize;
    for _ in 0..count {
        let id = number(tokens.next(), "element tag")?;
        let code: u32 = number(tokens.next(), "element type")?;
        let cell_kind = kind(code)?;
        let tag_count: usize = number(tokens.next(), "element tag count")?;
        if tag_count > source.len() / 2 {
            return Err(err("element tags exceed available input"));
        }
        let mut meaningful_tag = false;
        let mut physical = None;
        for index in 0..tag_count {
            let tag: i64 = number(tokens.next(), "element metadata tag")?;
            if index == 0 {
                if tag < 0 {
                    return Err(err("physical tag must be nonnegative"));
                }
                if tag > 0 {
                    physical = Some(u64::try_from(tag).map_err(|_| err("invalid physical tag"))?);
                }
            } else {
                meaningful_tag |= tag != 0;
            }
        }
        if meaningful_tag {
            tagged_elements += 1;
        }
        let mut nodes = Vec::with_capacity(cell_kind.node_count());
        for _ in 0..cell_kind.node_count() {
            nodes.push(number(tokens.next(), "element node")?);
        }
        parsed.push(Element22 {
            id,
            code,
            nodes,
            physical,
        });
    }
    if tokens.next().is_some() {
        return Err(err("extra MSH 2.2 element tokens"));
    }
    let names = physical_names(source)?;
    let mut selected_cells: BTreeMap<(u8, u64), Vec<u64>> = BTreeMap::new();
    for element in &parsed {
        if let Some(tag) = element.physical {
            let dimension = kind(element.code)?.dimension();
            selected_cells
                .entry((dimension, tag))
                .or_default()
                .push(element.id);
        }
    }

    // The 4.1 parser already owns node-ID resolution and complete data-block
    // validation. Translate only the flat geometry grammar, preserving tags.
    let mut text = String::from("$MeshFormat\n4.1 0 8\n$EndMeshFormat\n$Nodes\n");
    let min = points.iter().map(|(id, _)| *id).min().unwrap_or(0);
    let max = points.iter().map(|(id, _)| *id).max().unwrap_or(0);
    if points.is_empty() {
        text.push_str("0 0 0 0\n");
    } else {
        writeln!(
            text,
            "1 {} {min} {max}\n3 1 0 {}",
            points.len(),
            points.len()
        )
        .expect("String write");
        for (id, _) in &points {
            writeln!(text, "{id}").expect("String write");
        }
        for (_, xyz) in &points {
            writeln!(text, "{} {} {}", xyz[0], xyz[1], xyz[2]).expect("String write");
        }
    }
    text.push_str("$EndNodes\n$Elements\n");
    let mut groups: Vec<Vec<Element22>> = Vec::new();
    for element in parsed {
        if groups
            .last()
            .is_none_or(|group| group[0].code != element.code)
        {
            groups.push(Vec::new());
        }
        groups.last_mut().expect("group just created").push(element);
    }
    let element_count: usize = groups.iter().map(Vec::len).sum();
    let min = groups
        .iter()
        .flat_map(|group| group.iter().map(|item| item.id))
        .min()
        .unwrap_or(0);
    let max = groups
        .iter()
        .flat_map(|group| group.iter().map(|item| item.id))
        .max()
        .unwrap_or(0);
    writeln!(text, "{} {element_count} {min} {max}", groups.len()).expect("String write");
    for group in groups {
        let cell_kind = kind(group[0].code)?;
        writeln!(
            text,
            "{} 1 {} {}",
            code(cell_kind).1,
            group[0].code,
            group.len()
        )
        .expect("String write");
        for element in group {
            write!(text, "{}", element.id).expect("String write");
            for node in element.nodes {
                write!(text, " {node}").expect("String write");
            }
            text.push('\n');
        }
    }
    text.push_str("$EndElements\n");
    for name in ["NodeData", "ElementData"] {
        for body in data_sections(source, name)? {
            writeln!(text, "${name}\n{}\n$End{name}", body.trim()).expect("String write");
        }
    }
    let (mut dataset, _) = read_41(&text)?;
    let generated_group_names = attach_cell_sets(&mut dataset.mesh, &names, selected_cells);
    dataset.validate()?;
    Ok(Projection {
        dataset,
        tagged_elements,
        generated_group_names,
    })
}

/// Read ASCII Gmsh MSH 4.1 or 2.2 and report unrepresented 2.2 tags.
///
/// # Errors
///
/// Returns an error for malformed, binary, unsupported, or incomplete data.
pub fn read_projection(source: &str) -> Result<Projection> {
    let normalized = source.replace("\r\n", "\n");
    let header = section(&normalized, "MeshFormat")?.ok_or_else(|| err("missing MeshFormat"))?;
    let mut tokens = header.split_whitespace();
    let version = tokens.next().ok_or_else(|| err("missing MSH version"))?;
    if tokens.next() != Some("0") || tokens.next() != Some("8") || tokens.next().is_some() {
        return Err(err("only ASCII MSH with 8-byte data size is supported"));
    }
    match version {
        "4.1" => {
            let (dataset, generated_group_names) = read_41(&normalized)?;
            Ok(Projection {
                dataset,
                tagged_elements: 0,
                generated_group_names,
            })
        }
        "2.2" => read_22(&normalized),
        _ => Err(err(format!("unsupported MSH version {version}"))),
    }
}

/// Read ASCII Gmsh MSH 4.1 or 2.2 into the common dataset.
///
/// # Errors
///
/// Returns an error for malformed, binary, unsupported, or incomplete data.
#[cfg(test)]
pub fn read(source: &str) -> Result<Dataset> {
    Ok(read_projection(source)?.dataset)
}