Skip to main content

teamcraft/simmesh/
frd.rs

1//! CalculiX results (`.frd`, ASCII): nodes, elements and nodal result blocks.
2
3use std::collections::HashMap;
4
5use anyhow::{bail, Context, Result};
6
7use crate::mesh::{Field, Mesh, MAX_FIELDS};
8
9use super::*;
10
11/// Nodes of the FRD element types 1–12 and their VTK cell type.
12pub(super) fn frd_element(ty: usize) -> Option<(usize, u8)> {
13    Some(match ty {
14        1 => (8, HEXAHEDRON),
15        2 => (6, WEDGE),
16        3 => (4, TETRA),
17        4 => (20, QUADRATIC_HEXAHEDRON),
18        5 => (15, QUADRATIC_WEDGE),
19        6 => (10, QUADRATIC_TETRA),
20        7 => (3, TRIANGLE),
21        8 => (6, QUADRATIC_TRIANGLE),
22        9 => (4, QUAD),
23        10 => (8, QUADRATIC_QUAD),
24        11 => (2, LINE),
25        12 => (3, QUADRATIC_EDGE),
26        _ => return None,
27    })
28}
29
30/// Fixed-width record of an FRD data line: `' -1'`, the id in `width`
31/// columns, then numbers in 12 columns (E12.5). Falls back to whitespace
32/// separated tokens for files written by other programs.
33pub(super) fn frd_record(
34    line: &str,
35    width: usize,
36    values: usize,
37    out: &mut Vec<f64>,
38) -> Option<i64> {
39    let start = out.len();
40    let fixed = (|| {
41        let id = line.get(3..3 + width)?.trim();
42        let id = if id.is_empty() { 0 } else { id.parse().ok()? };
43        for k in 0..values {
44            let at = 3 + width + 12 * k;
45            let s = line.get(at..(at + 12).min(line.len()))?.trim();
46            out.push(s.parse().ok()?);
47        }
48        Some(id)
49    })();
50    if fixed.is_some() {
51        return fixed;
52    }
53    out.truncate(start);
54    let mut tokens = line.split_ascii_whitespace();
55    let tag = tokens.next()?;
56    let id = if tag == "-1" {
57        tokens.next()?.parse().ok()?
58    } else {
59        0
60    };
61    for _ in 0..values {
62        out.push(tokens.next()?.parse().ok()?);
63    }
64    Some(id)
65}
66
67/// Integers of an FRD element line, `width` columns each from column 3.
68pub(super) fn frd_ints(line: &str, width: usize, out: &mut Vec<i64>) -> Option<()> {
69    let mut at = 3;
70    let mut any = false;
71    while at < line.len() {
72        let s = line.get(at..(at + width).min(line.len()))?.trim();
73        if !s.is_empty() {
74            out.push(s.parse().ok()?);
75            any = true;
76        }
77        at += width;
78    }
79    any.then_some(())
80}
81
82/// Priority of FRD result blocks in the field selector (stress first).
83pub(super) fn frd_priority(name: &str) -> u8 {
84    if name.starts_with("STRESS") {
85        0
86    } else if name.starts_with("DISP") {
87        1
88    } else {
89        2
90    }
91}
92
93/// CalculiX result files (`.frd`, ASCII). Nodes, elements and every nodal
94/// result block are read; a block that occurs several times (one per
95/// increment) shows its last increment.
96pub fn frd(data: &[u8]) -> Result<Mesh> {
97    let text = String::from_utf8_lossy(data);
98    let mut lines = text.lines().map(|l| l.trim_end_matches('\r'));
99    let mut grid = Grid::default();
100    let mut ids: HashMap<i64, u32> = HashMap::new();
101    let mut width = 10;
102    // (block name, field) in order of appearance
103    let mut fields: Vec<(String, Field)> = Vec::new();
104    let mut seen_nodes = false;
105    let mut buf: Vec<f64> = Vec::new();
106    let mut ints: Vec<i64> = Vec::new();
107    while let Some(line) = lines.next() {
108        let head = line.trim_start();
109        if head.starts_with("2C") || head.starts_with("3C") {
110            let tokens: Vec<&str> = head.split_ascii_whitespace().collect();
111            let format = if tokens.len() >= 3 {
112                tokens.last().and_then(|t| t.parse::<u8>().ok())
113            } else {
114                None
115            };
116            match format {
117                Some(2) => bail!("Binary FRD files have no preview – write the results as ASCII (*NODE FILE, OUTPUT=ASCII)."),
118                Some(0) => width = 5,
119                Some(_) => width = 10,
120                None => {}
121            }
122            if head.starts_with("2C") {
123                if seen_nodes {
124                    bail!("the FRD file has two node blocks");
125                }
126                seen_nodes = true;
127                for line in lines.by_ref() {
128                    let tag = line.trim_start();
129                    if tag.starts_with("-3") {
130                        break;
131                    }
132                    if !tag.starts_with("-1") {
133                        bail!("unexpected line in the FRD node block");
134                    }
135                    buf.clear();
136                    let id =
137                        frd_record(line, width, 3, &mut buf).context("invalid FRD node line")?;
138                    if ids.len() >= MAX_NODES {
139                        bail!("The result has more than {MAX_NODES} nodes – too many for the preview.");
140                    }
141                    let next = grid.point_count() as u32;
142                    match ids.entry(id) {
143                        std::collections::hash_map::Entry::Occupied(o) => {
144                            // a node written twice: the later coordinates win
145                            let at = *o.get() as usize * 3;
146                            let point = grid.points.get_mut(at..at + 3).unwrap_or_default();
147                            for (p, &c) in point.iter_mut().zip(&buf) {
148                                *p = c as f32;
149                            }
150                        }
151                        std::collections::hash_map::Entry::Vacant(v) => {
152                            v.insert(next);
153                            grid.points.extend(buf.iter().map(|&c| c as f32));
154                        }
155                    }
156                }
157            } else {
158                // elements: ' -1' id type group material, then ' -2' node lines
159                let mut current: Option<(usize, u8)> = None;
160                ints.clear();
161                let finish = |grid: &mut Grid,
162                              current: Option<(usize, u8)>,
163                              ints: &mut Vec<i64>|
164                 -> Result<()> {
165                    if let Some((nodes, ty)) = current {
166                        let element = ints
167                            .get(..nodes)
168                            .context("an FRD element has too few nodes")?;
169                        grid.check_cells(1, nodes)?;
170                        for id in element {
171                            let p = ids
172                                .get(id)
173                                .context("an FRD element refers to a missing node")?;
174                            grid.conn.push(*p);
175                        }
176                        grid.ends.push(grid.conn.len());
177                        grid.types.push(ty);
178                    }
179                    ints.clear();
180                    Ok(())
181                };
182                for line in lines.by_ref() {
183                    let tag = line.trim_start();
184                    if tag.starts_with("-3") {
185                        break;
186                    } else if tag.starts_with("-1") {
187                        finish(&mut grid, current.take(), &mut ints)?;
188                        // id (width) then type (I5): read the type column directly
189                        let ty = line
190                            .get(3 + width..3 + width + 5)
191                            .and_then(|s| s.trim().parse::<usize>().ok())
192                            .or_else(|| {
193                                line.split_ascii_whitespace()
194                                    .nth(2)
195                                    .and_then(|s| s.parse().ok())
196                            })
197                            .context("invalid FRD element line")?;
198                        // unknown types (springs, …) are read and not drawn
199                        current = Some(frd_element(ty).unwrap_or((0, 0)));
200                    } else if tag.starts_with("-2") {
201                        if frd_ints(line, width, &mut ints).is_none() {
202                            bail!("invalid FRD element node line");
203                        }
204                        if ints.len() > 64 {
205                            bail!("an FRD element has too many nodes");
206                        }
207                    } else {
208                        bail!("unexpected line in the FRD element block");
209                    }
210                }
211                finish(&mut grid, current.take(), &mut ints)?;
212            }
213        } else if head.starts_with("-4") {
214            // result block: -4 NAME ncomps irtype, -5 component lines, data
215            let tokens: Vec<&str> = head.split_ascii_whitespace().collect();
216            let name = tokens.get(1).copied().unwrap_or("RESULT").to_string();
217            let comps = count(tokens.get(2).copied(), "FRD component count")?;
218            if comps > MAX_COMPONENTS {
219                bail!("an FRD result block has too many components");
220            }
221            let mut stored = 0usize;
222            let mut tensor = true;
223            for _ in 0..comps {
224                let c = lines.next().context("the FRD file is truncated")?;
225                let t: Vec<&str> = c.split_ascii_whitespace().collect();
226                if t.first() != Some(&"-5") {
227                    bail!("component line (-5) expected in the FRD result block");
228                }
229                // the sixth number (iexist) is 1 for components the post-processor computes
230                let computed = t.get(6).is_some_and(|s| s.starts_with('1'));
231                if !computed {
232                    stored += 1;
233                    tensor &= t.get(3) == Some(&"4");
234                }
235            }
236            let n = grid.point_count();
237            // At most MAX_FIELDS blocks are kept (a later increment replaces
238            // its block, a higher-priority block replaces a lower one); the
239            // others are read and not stored, so a file with thousands of
240            // differently named blocks cannot allocate a value array each.
241            let slot = match fields.iter().position(|(k, _)| *k == name) {
242                Some(i) => Some(i),
243                None if fields.len() < MAX_FIELDS => Some(fields.len()),
244                None => fields
245                    .iter()
246                    .rposition(|(k, _)| frd_priority(k) > frd_priority(&name)),
247            };
248            let mut values = if slot.is_some() {
249                vec![f32::NAN; n]
250            } else {
251                Vec::new()
252            };
253            let reduce = if stored == 6 && tensor {
254                Reduce::VonMises6
255            } else if stored <= 1 {
256                Reduce::Value
257            } else {
258                Reduce::Magnitude
259            };
260            let mut node: Option<u32> = None;
261            let flush = |node: Option<u32>, buf: &mut Vec<f64>, values: &mut Vec<f32>| {
262                if let Some(v) = node.and_then(|p| values.get_mut(p as usize)) {
263                    if let Some(components) = buf.get(..stored).filter(|_| stored > 0) {
264                        *v = reduce.apply(components);
265                    }
266                }
267                buf.clear();
268            };
269            buf.clear();
270            for line in lines.by_ref() {
271                let tag = line.trim_start();
272                if tag.starts_with("-3") {
273                    break;
274                } else if tag.starts_with("-1") {
275                    flush(node.take(), &mut buf, &mut values);
276                    let id = frd_record(line, width, stored.min(6), &mut buf)
277                        .context("invalid FRD result line")?;
278                    node = ids.get(&id).copied();
279                } else if tag.starts_with("-2") {
280                    let more = stored.saturating_sub(buf.len()).min(6);
281                    frd_record(line, width, more, &mut buf).context("invalid FRD result line")?;
282                } else {
283                    bail!("unexpected line in the FRD result block");
284                }
285            }
286            flush(node.take(), &mut buf, &mut values);
287            if let (Some(i), true) = (slot, stored > 0 && n > 0) {
288                let field = Field {
289                    name: reduce.label(&name),
290                    unit: String::new(),
291                    values,
292                };
293                if let Some((_, slot)) = fields.get_mut(i).filter(|(k, _)| *k == name) {
294                    *slot = field;
295                } else {
296                    if i < fields.len() {
297                        fields.remove(i); // a lower-priority block makes room
298                    }
299                    fields.push((name, field));
300                }
301            }
302        } else if head.starts_with("9999") {
303            break;
304        }
305    }
306    if !seen_nodes {
307        bail!("Not a CalculiX result file (no node block).");
308    }
309    fields.sort_by_key(|(name, _)| frd_priority(name));
310    for (_, f) in fields {
311        Grid::add_field(&mut grid.point_fields, f);
312    }
313    grid.into_mesh()
314}
315
316#[cfg(test)]
317mod unit_tests {
318    use super::*;
319
320    #[test]
321    fn calculix_element_types() {
322        let expected = [
323            (1, 8, HEXAHEDRON),
324            (2, 6, WEDGE),
325            (3, 4, TETRA),
326            (4, 20, QUADRATIC_HEXAHEDRON),
327            (5, 15, QUADRATIC_WEDGE),
328            (6, 10, QUADRATIC_TETRA),
329            (7, 3, TRIANGLE),
330            (8, 6, QUADRATIC_TRIANGLE),
331            (9, 4, QUAD),
332            (10, 8, QUADRATIC_QUAD),
333            (11, 2, LINE),
334            (12, 3, QUADRATIC_EDGE),
335        ];
336        for (ty, nodes, vtk) in expected {
337            assert_eq!(frd_element(ty), Some((nodes, vtk)), "type {ty}");
338        }
339        assert_eq!(frd_element(0), None);
340        assert_eq!(frd_element(13), None);
341    }
342
343    #[test]
344    fn fixed_width_records() {
345        // " -1" + id (10 wide) + values (12 wide each), as CalculiX writes them
346        let line = format!(
347            " -1{:>10}{:>12}{:>12}{:>12}",
348            7, "1.50000E+00", "-2.0E-01", "3.00000E+00"
349        );
350        let mut out = Vec::new();
351        assert_eq!(frd_record(&line, 10, 3, &mut out), Some(7));
352        assert_eq!(out, [1.5, -0.2, 3.0]);
353        // values that touch each other only work by width
354        let packed = " -1         1-1.00000E+00-2.00000E+00";
355        let mut out = Vec::new();
356        assert_eq!(frd_record(packed, 10, 2, &mut out), Some(1));
357        assert_eq!(out, [-1.0, -2.0]);
358        // a free-format line is read by tokens
359        let mut out = vec![9.0];
360        assert_eq!(frd_record(" -1 12 4.5 5.5", 5, 2, &mut out), Some(12));
361        assert_eq!(out, [9.0, 4.5, 5.5], "appended after what was there");
362        // a continuation line without an id
363        let mut out = Vec::new();
364        assert_eq!(frd_record(" -2 6.5", 5, 1, &mut out), Some(0));
365        assert_eq!(out, [6.5]);
366        // too few values: the line is refused (every caller then stops with
367        // an error, so values read before the gap are never used)
368        let mut out = Vec::new();
369        assert_eq!(frd_record(" -1 3 4.5", 5, 3, &mut out), None);
370        assert_eq!(frd_record("", 10, 1, &mut Vec::new()), None);
371    }
372
373    #[test]
374    fn fixed_width_integers() {
375        let mut out = Vec::new();
376        assert_eq!(
377            frd_ints(" -2         1         2         3", 10, &mut out),
378            Some(())
379        );
380        assert_eq!(out, [1, 2, 3]);
381        let mut out = Vec::new();
382        assert_eq!(frd_ints(" -2    ", 10, &mut out), None, "no numbers");
383        assert_eq!(frd_ints(" -2      abc", 10, &mut out), None, "not a number");
384    }
385
386    #[test]
387    fn stress_first_then_displacements() {
388        assert_eq!(frd_priority("STRESS"), 0);
389        assert_eq!(frd_priority("STRESSR"), 0);
390        assert_eq!(frd_priority("DISP"), 1);
391        assert_eq!(frd_priority("NDTEMP"), 2);
392        assert_eq!(frd_priority("stress"), 2, "names are upper case in FRD");
393    }
394}