enok/src/inndata/ifcgeom.rs

//! Geometri fra IFC: arealer og flatenormaler fra trekantnett, brep og ekstruderte profiler.
//!
//! BIM-verktøy eksporterer sjelden mengder for alle bygningsdeler, og IFC-visningen
//! ReferenceView bruker tessellert geometri. Denne modulen regner derfor arealet av en
//! bygningsdel direkte av geometrien: trekantene grupperes etter retning, og den største
//! gruppen av parallelle flater er delens hovedflate. Summen deles på to fordi et lukket
//! volum har en forside og en bakside.
//!
//! Metoden gir riktig areal også når vinduer og dører er skåret ut av veggen, siden hullene
//! rett og slett mangler trekanter.

use crate::inndata::step::{Arg, Fil};

/// Et trekantnett i bygningsdelens lokale koordinater.
#[derive(Default, Debug)]
pub struct Nett {
    pub punkter: Vec<[f64; 3]>,
    pub trekanter: Vec<[usize; 3]>,
}

/// Hovedflaten til en bygningsdel.
#[derive(Debug, Clone, Copy)]
pub struct Hovedflate {
    /// Areal av én side, m².
    pub areal: f64,
    /// Enhetsnormal i lokale koordinater.
    pub normal: [f64; 3],
    /// Utstrekning i flatens plan: (u_min, u_maks, v_min, v_maks), m.
    plan: [f64; 4],
}

impl Nett {
    pub fn tom(&self) -> bool {
        self.trekanter.is_empty()
    }

    fn trekant(&self, t: [usize; 3]) -> Option<([f64; 3], f64)> {
        let (a, b, c) = (
            *self.punkter.get(t[0])?,
            *self.punkter.get(t[1])?,
            *self.punkter.get(t[2])?,
        );
        let u = [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
        let v = [c[0] - a[0], c[1] - a[1], c[2] - a[2]];
        let n = [
            u[1] * v[2] - u[2] * v[1],
            u[2] * v[0] - u[0] * v[2],
            u[0] * v[1] - u[1] * v[0],
        ];
        let lengde = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
        if lengde < 1e-12 {
            return None;
        }
        Some(([n[0] / lengde, n[1] / lengde, n[2] / lengde], 0.5 * lengde))
    }

    /// Samlet overflate, m². Tar med begge sider av et lukket volum.
    pub fn overflate(&self, skala: f64) -> f64 {
        let s2 = skala * skala;
        self.trekanter
            .iter()
            .filter_map(|t| self.trekant(*t))
            .map(|(_, a)| a * s2)
            .sum()
    }

    /// Finner den største gruppen av parallelle trekanter, som er bygningsdelens hovedflate.
    ///
    /// `skala` gjør lengdene om til meter. Trekanter regnes som parallelle når normalene
    /// er innenfor omtrent 18 grader, målt uten hensyn til retning, slik at for- og bakside
    /// havner i samme gruppe.
    pub fn hovedflate(&self, skala: f64) -> Option<Hovedflate> {
        let s2 = skala * skala;
        let mut flater: Vec<([f64; 3], f64)> = self
            .trekanter
            .iter()
            .filter_map(|t| self.trekant(*t))
            .map(|(n, a)| (n, a * s2))
            .collect();
        if flater.is_empty() {
            return None;
        }
        flater.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));

        let mut grupper: Vec<([f64; 3], f64)> = Vec::new();
        for (n, a) in &flater {
            let treff = grupper
                .iter_mut()
                .find(|(g, _)| (g[0] * n[0] + g[1] * n[1] + g[2] * n[2]).abs() > 0.95);
            match treff {
                Some((_, sum)) => *sum += a,
                None => grupper.push((*n, *a)),
            }
        }
        let (normal, sum) = grupper
            .into_iter()
            .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))?;
        Some(Hovedflate {
            areal: sum / 2.0,
            normal,
            plan: self.utstrekning(normal, skala),
        })
    }

    /// Punktenes utstrekning i planet vinkelrett på `normal`, m.
    fn utstrekning(&self, normal: [f64; 3], skala: f64) -> [f64; 4] {
        let (u, v) = planakser(normal);
        let mut ut = [
            f64::INFINITY,
            f64::NEG_INFINITY,
            f64::INFINITY,
            f64::NEG_INFINITY,
        ];
        for p in &self.punkter {
            let a = (p[0] * u[0] + p[1] * u[1] + p[2] * u[2]) * skala;
            let b = (p[0] * v[0] + p[1] * v[1] + p[2] * v[2]) * skala;
            ut[0] = ut[0].min(a);
            ut[1] = ut[1].max(a);
            ut[2] = ut[2].min(b);
            ut[3] = ut[3].max(b);
        }
        ut
    }
}

/// To akser som spenner planet vinkelrett på `n`.
fn planakser(n: [f64; 3]) -> ([f64; 3], [f64; 3]) {
    let hjelp = if n[2].abs() < 0.9 {
        [0.0, 0.0, 1.0]
    } else {
        [1.0, 0.0, 0.0]
    };
    let u = enhet([
        n[1] * hjelp[2] - n[2] * hjelp[1],
        n[2] * hjelp[0] - n[0] * hjelp[2],
        n[0] * hjelp[1] - n[1] * hjelp[0],
    ])
    .unwrap_or([1.0, 0.0, 0.0]);
    let v = [
        n[1] * u[2] - n[2] * u[1],
        n[2] * u[0] - n[0] * u[2],
        n[0] * u[1] - n[1] * u[0],
    ];
    (u, v)
}

/// Slår sammen hovedflatene til en bygningsdels geometrideler til ett areal.
///
/// BIM-verktøy eksporterer ofte hvert materialsjikt i en vegg eller et dekke som sin egen
/// geometridel. Sjiktene ligger parallelt og dekker den samme flaten, så de skal telles én
/// gang. Deler som ligger i samme plan og overlapper hverandre regnes derfor som én flate,
/// der den største av dem gjelder, mens deler som dekker forskjellige områder legges sammen.
pub fn samlet_hovedflate(deler: &[Nett], skala: f64) -> Option<Hovedflate> {
    let mut flater: Vec<Hovedflate> = deler.iter().filter_map(|n| n.hovedflate(skala)).collect();
    if flater.is_empty() {
        return None;
    }
    flater.sort_by(|a, b| {
        b.areal
            .partial_cmp(&a.areal)
            .unwrap_or(std::cmp::Ordering::Equal)
    });
    let mut grupper: Vec<Hovedflate> = Vec::new();
    for f in flater {
        let samme = grupper.iter().any(|g| {
            let parallell =
                (g.normal[0] * f.normal[0] + g.normal[1] * f.normal[1] + g.normal[2] * f.normal[2])
                    .abs()
                    > 0.95;
            parallell && overlapper(&g.plan, &f.plan)
        });
        if !samme {
            grupper.push(f);
        }
    }
    let areal = grupper.iter().map(|g| g.areal).sum();
    Some(Hovedflate {
        areal,
        ..grupper[0]
    })
}

/// Sant når to utstrekninger dekker mer enn halvparten av den minste i begge retninger.
fn overlapper(a: &[f64; 4], b: &[f64; 4]) -> bool {
    let dekning = |a0: f64, a1: f64, b0: f64, b1: f64| {
        let felles = (a1.min(b1) - a0.max(b0)).max(0.0);
        let minste = (a1 - a0).min(b1 - b0);
        if minste <= 1e-9 {
            felles > 0.0
        } else {
            felles / minste > 0.5
        }
    };
    dekning(a[0], a[1], b[0], b[1]) && dekning(a[2], a[3], b[2], b[3])
}

/// Henter geometrien til et produkt via `Representation`, med én oppføring per geometridel.
pub fn netter_for_produkt(fil: &Fil, representasjon: Option<u64>) -> Vec<Nett> {
    let mut deler = Vec::new();
    let Some(r) = representasjon else {
        return deler;
    };
    let Some(form) = fil.hent(r) else {
        return deler;
    };
    // IfcProductDefinitionShape.Representations
    for rep in form.arg(2).map(|a| a.som_liste()).unwrap_or(&[]) {
        let Some(sr) = rep.som_ref().and_then(|x| fil.hent(x)) else {
            continue;
        };
        if sr.typ != "IFCSHAPEREPRESENTATION" {
            continue;
        }
        // Bare kroppen beskriver volumet; akse- og fotavtrykksrepresentasjoner hoppes over.
        if sr.arg(1).and_then(|a| a.som_tekst()) != Some("Body") {
            continue;
        }
        for post in sr.arg(3).map(|a| a.som_liste()).unwrap_or(&[]) {
            if let Some(n) = post.som_ref() {
                let mut nett = Nett::default();
                les_element(fil, n, &mut nett, 0);
                if !nett.tom() {
                    deler.push(nett);
                }
            }
        }
    }
    deler
}

fn les_element(fil: &Fil, r: u64, nett: &mut Nett, dybde: u32) {
    if dybde > 8 {
        return;
    }
    let Some(e) = fil.hent(r) else {
        return;
    };
    match e.typ.as_str() {
        "IFCTRIANGULATEDFACESET" => les_trekantnett(fil, r, nett),
        "IFCPOLYGONALFACESET" => les_polygonnett(fil, r, nett),
        "IFCFACETEDBREP" | "IFCADVANCEDBREP" => {
            if let Some(skall) = e.arg(0).and_then(|a| a.som_ref()) {
                les_skall(fil, skall, nett);
            }
        }
        "IFCEXTRUDEDAREASOLID" => les_ekstrudering(fil, r, nett),
        "IFCMAPPEDITEM" => {
            // Kartlagt geometri: følg kilden. Transformasjonen ses bort fra, siden
            // arealet er uavhengig av flytting og rotasjon.
            if let Some(kilde) = e.arg(0).and_then(|a| a.som_ref()).and_then(|k| fil.hent(k)) {
                if let Some(rep) = kilde
                    .arg(1)
                    .and_then(|a| a.som_ref())
                    .and_then(|x| fil.hent(x))
                {
                    for post in rep.arg(3).map(|a| a.som_liste()).unwrap_or(&[]) {
                        if let Some(n) = post.som_ref() {
                            les_element(fil, n, nett, dybde + 1);
                        }
                    }
                }
            }
        }
        "IFCBOOLEANRESULT" | "IFCBOOLEANCLIPPINGRESULT" => {
            // Bruk førsteoperanden; utsparinger gir uansett litt for stort areal.
            if let Some(f) = e.arg(1).and_then(|a| a.som_ref()) {
                les_element(fil, f, nett, dybde + 1);
            }
        }
        _ => {}
    }
}

fn punktliste(fil: &Fil, r: u64) -> Vec<[f64; 3]> {
    let Some(e) = fil.hent(r) else {
        return Vec::new();
    };
    e.arg(0)
        .map(|a| a.som_liste())
        .unwrap_or(&[])
        .iter()
        .map(|p| {
            let k = p.som_liste();
            let g = |i: usize| k.get(i).and_then(|x| x.som_tall()).unwrap_or(0.0);
            [g(0), g(1), g(2)]
        })
        .collect()
}

fn les_trekantnett(fil: &Fil, r: u64, nett: &mut Nett) {
    let Some(e) = fil.hent(r) else { return };
    let Some(liste) = e.arg(0).and_then(|a| a.som_ref()) else {
        return;
    };
    let start = nett.punkter.len();
    nett.punkter.extend(punktliste(fil, liste));
    // CoordIndex er 1-basert.
    for t in e.arg(3).map(|a| a.som_liste()).unwrap_or(&[]) {
        let k = t.som_liste();
        if k.len() < 3 {
            continue;
        }
        let g = |i: usize| k[i].som_tall().unwrap_or(1.0) as usize;
        nett.trekanter.push([
            start + g(0).saturating_sub(1),
            start + g(1).saturating_sub(1),
            start + g(2).saturating_sub(1),
        ]);
    }
}

fn les_polygonnett(fil: &Fil, r: u64, nett: &mut Nett) {
    let Some(e) = fil.hent(r) else { return };
    let Some(liste) = e.arg(0).and_then(|a| a.som_ref()) else {
        return;
    };
    let start = nett.punkter.len();
    nett.punkter.extend(punktliste(fil, liste));
    for f in e.arg(2).map(|a| a.som_liste()).unwrap_or(&[]) {
        let Some(fe) = f.som_ref().and_then(|x| fil.hent(x)) else {
            continue;
        };
        let indekser: Vec<usize> = fe
            .arg(0)
            .map(|a| a.som_liste())
            .unwrap_or(&[])
            .iter()
            .filter_map(|x| x.som_tall())
            .map(|v| start + (v as usize).saturating_sub(1))
            .collect();
        vifte(&indekser, nett);
    }
}

fn les_skall(fil: &Fil, r: u64, nett: &mut Nett) {
    let Some(skall) = fil.hent(r) else { return };
    for f in skall.arg(0).map(|a| a.som_liste()).unwrap_or(&[]) {
        let Some(flate) = f.som_ref().and_then(|x| fil.hent(x)) else {
            continue;
        };
        for b in flate.arg(0).map(|a| a.som_liste()).unwrap_or(&[]) {
            let Some(grense) = b.som_ref().and_then(|x| fil.hent(x)) else {
                continue;
            };
            // Bare ytre grense; hull i flaten gir litt for stort areal.
            if grense.typ != "IFCFACEOUTERBOUND" {
                continue;
            }
            let Some(lokke) = grense
                .arg(0)
                .and_then(|a| a.som_ref())
                .and_then(|x| fil.hent(x))
            else {
                continue;
            };
            let mut indekser = Vec::new();
            for p in lokke.arg(0).map(|a| a.som_liste()).unwrap_or(&[]) {
                if let Some(pr) = p.som_ref() {
                    if let Some(pe) = fil.hent(pr) {
                        let k = pe.arg(0).map(|a| a.som_liste()).unwrap_or(&[]);
                        let g = |i: usize| k.get(i).and_then(|x| x.som_tall()).unwrap_or(0.0);
                        indekser.push(nett.punkter.len());
                        nett.punkter.push([g(0), g(1), g(2)]);
                    }
                }
            }
            vifte(&indekser, nett);
        }
    }
}

/// Deler en polygonløkke i trekanter fra første hjørne.
fn vifte(indekser: &[usize], nett: &mut Nett) {
    for i in 1..indekser.len().saturating_sub(1) {
        nett.trekanter
            .push([indekser[0], indekser[i], indekser[i + 1]]);
    }
}

fn les_ekstrudering(fil: &Fil, r: u64, nett: &mut Nett) {
    let Some(e) = fil.hent(r) else { return };
    let Some(profil) = e.arg(0).and_then(|a| a.som_ref()).and_then(|x| fil.hent(x)) else {
        return;
    };
    let dybde = e.arg(3).and_then(|a| a.som_tall()).unwrap_or(0.0);
    let retning = e
        .arg(2)
        .and_then(|a| a.som_ref())
        .and_then(|x| fil.hent(x))
        .map(|d| {
            let k = d.arg(0).map(|a| a.som_liste()).unwrap_or(&[]);
            let g = |i: usize| k.get(i).and_then(|x| x.som_tall()).unwrap_or(0.0);
            [g(0), g(1), g(2)]
        })
        .unwrap_or([0.0, 0.0, 1.0]);

    // Profilet ligger i xy-planet. Hent ytre kurve som punkter.
    let punkter = profilpunkter(fil, profil);
    if punkter.len() < 3 {
        return;
    }
    let start = nett.punkter.len();
    for p in &punkter {
        nett.punkter.push([p[0], p[1], 0.0]);
    }
    let topp = nett.punkter.len();
    for p in &punkter {
        nett.punkter.push([
            p[0] + retning[0] * dybde,
            p[1] + retning[1] * dybde,
            retning[2] * dybde,
        ]);
    }
    // Bunn og topp.
    let bunn: Vec<usize> = (start..topp).collect();
    let topp_i: Vec<usize> = (topp..nett.punkter.len()).collect();
    vifte(&bunn, nett);
    vifte(&topp_i, nett);
    // Sideflater.
    let n = punkter.len();
    for i in 0..n {
        let j = (i + 1) % n;
        nett.trekanter.push([start + i, start + j, topp + j]);
        nett.trekanter.push([start + i, topp + j, topp + i]);
    }
}

fn profilpunkter(fil: &Fil, profil: &crate::inndata::step::Entitet) -> Vec<[f64; 2]> {
    match profil.typ.as_str() {
        "IFCARBITRARYCLOSEDPROFILEDEF" | "IFCARBITRARYPROFILEDEFWITHVOIDS" => {
            let Some(kurve) = profil
                .arg(2)
                .and_then(|a| a.som_ref())
                .and_then(|x| fil.hent(x))
            else {
                return Vec::new();
            };
            if kurve.typ != "IFCPOLYLINE" {
                return Vec::new();
            }
            kurve
                .arg(0)
                .map(|a| a.som_liste())
                .unwrap_or(&[])
                .iter()
                .filter_map(|p| p.som_ref().and_then(|x| fil.hent(x)))
                .map(|pe| {
                    let k = pe.arg(0).map(|a| a.som_liste()).unwrap_or(&[]);
                    let g = |i: usize| k.get(i).and_then(|x| x.som_tall()).unwrap_or(0.0);
                    [g(0), g(1)]
                })
                .collect()
        }
        "IFCRECTANGLEPROFILEDEF" => {
            let x = profil.arg(3).and_then(|a| a.som_tall()).unwrap_or(0.0) / 2.0;
            let y = profil.arg(4).and_then(|a| a.som_tall()).unwrap_or(0.0) / 2.0;
            vec![[-x, -y], [x, -y], [x, y], [-x, y]]
        }
        _ => Vec::new(),
    }
}

/// Lengdeenheten i meter, hentet gjennom prosjektets enhetsoppsett.
///
/// En IFC-fil kan inneholde flere `IfcSIUnit` med lengdeenhet, for eksempel når et
/// underoppsett bruker meter mens modellen er i millimeter. Derfor må enheten slås opp
/// gjennom `IfcProject.UnitsInContext` og ikke ved å lete gjennom alle enheter.
pub fn lengdeenhet(fil: &Fil) -> f64 {
    let fra_oppsett = fil
        .av_type("IFCPROJECT")
        .next()
        .and_then(|(_, p)| p.arg(8).and_then(|a| a.som_ref()))
        .and_then(|r| fil.hent(r))
        .and_then(|oppsett| {
            oppsett
                .arg(0)
                .map(|a| a.som_liste())
                .unwrap_or(&[])
                .iter()
                .filter_map(|u| u.som_ref().and_then(|r| fil.hent(r)))
                .find(|e| {
                    e.typ == "IFCSIUNIT"
                        && e.arg(1).and_then(|a| a.som_tekst()) == Some("LENGTHUNIT")
                })
                .map(prefiksfaktor)
        });
    fra_oppsett.unwrap_or(1.0)
}

fn prefiksfaktor(e: &crate::inndata::step::Entitet) -> f64 {
    match e.arg(2).and_then(|a| a.som_tekst()).unwrap_or("") {
        "MILLI" => 0.001,
        "CENTI" => 0.01,
        "DECI" => 0.1,
        "KILO" => 1000.0,
        _ => 1.0,
    }
}

/// Enhetsvektor, eller `None` for en nullvektor.
pub fn enhet(v: [f64; 3]) -> Option<[f64; 3]> {
    let l = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt();
    (l > 1e-12).then(|| [v[0] / l, v[1] / l, v[2] / l])
}

/// Sant når et argument er en liste med minst ett element.
pub fn har_innhold(a: Option<&Arg>) -> bool {
    a.map(|x| !x.som_liste().is_empty()).unwrap_or(false)
}

#[cfg(test)]
mod tester {
    use super::*;

    /// En kasse på 4 × 3 × 0,2 m skal gi hovedflate 12 m² med normal langs z.
    #[test]
    fn hovedflate_av_kasse() {
        let mut n = Nett::default();
        let (x, y, z) = (4.0, 3.0, 0.2);
        let p = [
            [0.0, 0.0, 0.0],
            [x, 0.0, 0.0],
            [x, y, 0.0],
            [0.0, y, 0.0],
            [0.0, 0.0, z],
            [x, 0.0, z],
            [x, y, z],
            [0.0, y, z],
        ];
        n.punkter.extend_from_slice(&p);
        let flater = [
            [0, 1, 2],
            [0, 2, 3], // bunn
            [4, 6, 5],
            [4, 7, 6], // topp
            [0, 4, 5],
            [0, 5, 1],
            [1, 5, 6],
            [1, 6, 2],
            [2, 6, 7],
            [2, 7, 3],
            [3, 7, 4],
            [3, 4, 0],
        ];
        n.trekanter.extend_from_slice(&flater);
        let h = n.hovedflate(1.0).unwrap();
        assert!((h.areal - 12.0).abs() < 1e-9, "areal ble {}", h.areal);
        assert!(h.normal[2].abs() > 0.99, "normal ble {:?}", h.normal);
    }

    /// Et hull i flaten skal trekkes fra, siden hullet mangler trekanter.
    #[test]
    fn hull_reduserer_hovedflaten() {
        let mut n = Nett::default();
        // To flater på 4 × 3 med et hull på 1 × 1 er modellert som fire striper.
        for z in [0.0, 0.2] {
            let base = n.punkter.len();
            n.punkter.extend_from_slice(&[
                [0.0, 0.0, z],
                [4.0, 0.0, z],
                [4.0, 1.0, z],
                [0.0, 1.0, z],
            ]);
            n.trekanter.push([base, base + 1, base + 2]);
            n.trekanter.push([base, base + 2, base + 3]);
        }
        let h = n.hovedflate(1.0).unwrap();
        assert!((h.areal - 4.0).abs() < 1e-9, "areal ble {}", h.areal);
    }

    #[test]
    fn skala_gjor_om_millimeter_til_kvadratmeter() {
        let mut n = Nett::default();
        n.punkter.extend_from_slice(&[
            [0.0, 0.0, 0.0],
            [2000.0, 0.0, 0.0],
            [2000.0, 1000.0, 0.0],
            [0.0, 1000.0, 0.0],
        ]);
        n.trekanter.push([0, 1, 2]);
        n.trekanter.push([0, 2, 3]);
        // Én enkelt flate: hovedflaten deler på to, så 2 m² blir 1 m².
        let h = n.hovedflate(0.001).unwrap();
        assert!((h.areal - 1.0).abs() < 1e-9, "areal ble {}", h.areal);
    }
}