conspire 0.7.7

The Rust interface to conspire.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

use crate::geometry::ntree::node::slot::Slot;
#[cfg(test)]
mod test;

use crate::{
    geometry::{
        Coordinate,
        ntree::{
            Balancing, Orthotree, Pairing, Rescaling,
            node::{Kind, Node, cell::Cell},
        },
    },
    io::{
        invalid,
        read::{attribute, bits, data_arrays, encoding, find_data_array, floats, tag},
        unsupported,
    },
    math::{Quantity, Scalar},
};
use std::{array::from_fn, collections::VecDeque, fs::read_to_string, io::Result, path::Path};

pub(crate) trait ReadHtg<P>
where
    P: AsRef<Path>,
    Self: Sized,
{
    fn read_htg(input: P) -> Result<Self>;
}

impl<const D: usize, const L: usize, const M: usize, const N: usize, T, U, P> ReadHtg<P>
    for Orthotree<D, L, M, N, T, U>
where
    P: AsRef<Path>,
    T: Cell,
    U: Slot,
{
    fn read_htg(input: P) -> Result<Self> {
        let text = read_to_string(input)?;
        let header = tag(&text, "<VTKFile")?;
        if attribute(header, "type") != Some("HyperTreeGrid") {
            return Err(invalid("not a HyperTreeGrid".into()));
        }
        if matches!(attribute(header, "byte_order"), Some(order) if order != "LittleEndian") {
            return Err(unsupported("big-endian HTG is not supported"));
        }
        let encoding = encoding(header)?;
        let grid = tag(&text, "<HyperTreeGrid")?;
        if matches!(attribute(grid, "BranchFactor"), Some(factor) if factor != "2") {
            return Err(unsupported("only BranchFactor 2 is supported"));
        }
        let dimensions =
            attribute(grid, "Dimensions").ok_or_else(|| invalid("no Dimensions".into()))?;
        let split_axes = dimensions.split_whitespace().filter(|d| *d != "1").count();
        if split_axes != D {
            return Err(invalid(format!(
                "HTG has {split_axes} refined axes but Orthotree was asked for D={D}"
            )));
        }
        let arrays = data_arrays(&text)?;
        let axes = ["XCoordinates", "YCoordinates", "ZCoordinates"];
        let mut lo = [0.0; D];
        let mut hi = [0.0; D];
        for (axis, &name) in axes.iter().take(D).enumerate() {
            let coordinates = floats(&find_data_array(&arrays, Some(name))?, &encoding)?;
            lo[axis] = coordinates[0] as Scalar;
            hi[axis] = coordinates[1] as Scalar;
        }
        let levels: usize = attribute(tag(&text, "<Tree ")?, "NumberOfLevels")
            .and_then(|n| n.parse().ok())
            .ok_or_else(|| invalid("no NumberOfLevels".into()))?;
        let descriptor = bits(&find_data_array(&arrays, Some("Descriptor"))?, &encoding)?;
        let root_length = 1usize << (levels - 1);
        let cell = Quantity::new((hi[0] - lo[0]) / root_length as Scalar);
        let rescale = Rescaling {
            center: Coordinate::from(from_fn(|axis| 0.5 * (lo[axis] + hi[axis]))),
            cell,
            half: 0.5 * root_length as Scalar,
        };
        let zero = number(0)?;
        let mut tree = Orthotree {
            balanced: Balancing::None,
            nodes: vec![Node {
                corner: from_fn(|_| zero),
                length: number(root_length)?,
                facets: [None; M],
                kind: Kind::Leaf,
                value: None,
            }],
            paired: Pairing::None,
            rescale,
        };
        let mut queue = VecDeque::<usize>::from([0]);
        let mut bit = 0;
        while let Some(node) = queue.pop_front() {
            if bit >= descriptor.len() {
                break;
            }
            let refined = descriptor[bit] == 1;
            bit += 1;
            if refined {
                tree.subdivide(node).map_err(|e| invalid(e.into()))?;
                queue.extend(
                    tree.nodes[node]
                        .orthants()
                        .unwrap()
                        .iter()
                        .map(|node: &U| node.slot()),
                );
            }
        }
        Ok(tree)
    }
}

fn number<T: Cell>(value: usize) -> Result<T> {
    T::length(value).ok_or_else(|| invalid("tree coordinate does not fit in T".into()))
}