1use std::fs;
13use std::path::{Path, PathBuf};
14
15use ndarray::{Array1, Array2};
16use num_complex::Complex64;
17
18use crate::core::types::{
19 BoundaryCondition, Element, ElementProperty, ElementType, PhysicsParams,
20};
21
22#[derive(Debug, Clone)]
24pub struct NcInputConfig {
25 pub version: String,
27 pub description: String,
29 pub control_params_i: Vec<i32>,
31 pub control_params_ii: Vec<f64>,
33 pub frequency_curve: Vec<(f64, f64, f64)>,
35 pub main_params_i: MainParamsI,
37 pub main_params_ii: MainParamsII,
39 pub main_params_iii: Vec<i32>,
41 pub main_params_iv: MainParamsIV,
43 pub node_files: Vec<PathBuf>,
45 pub element_files: Vec<PathBuf>,
47 pub symmetry: Option<SymmetryConfig>,
49 pub boundary_conditions: Vec<BoundarySpec>,
51 pub plane_waves: Vec<PlaneWaveSource>,
53 pub point_sources: Vec<PointSource>,
55 pub base_dir: PathBuf,
57}
58
59#[derive(Debug, Clone, Default)]
61pub struct MainParamsI {
62 pub element_type: i32,
64 pub num_nodes: usize,
66 pub num_elements: usize,
68 pub num_object_files: i32,
70 pub num_eval_files: i32,
72 pub bc_type: i32,
74 pub solver_method: i32,
76 pub fmm_method: i32,
78 pub parallel: i32,
80}
81
82#[derive(Debug, Clone, Default)]
84pub struct MainParamsII {
85 pub preconditioner: i32,
87 pub iterative_solver: i32,
89 pub reserved1: i32,
91 pub reserved2: f64,
93 pub output_level: i32,
95 pub reserved3: i32,
97 pub reserved4: i32,
99}
100
101#[derive(Debug, Clone)]
103pub struct MainParamsIV {
104 pub speed_of_sound: f64,
106 pub density: f64,
108 pub reference_pressure: f64,
110 pub reserved: Vec<f64>,
112}
113
114impl Default for MainParamsIV {
115 fn default() -> Self {
116 Self {
117 speed_of_sound: 343.0,
118 density: 1.21,
119 reference_pressure: 1.0,
120 reserved: vec![0.0; 4],
121 }
122 }
123}
124
125#[derive(Debug, Clone)]
127pub struct SymmetryConfig {
128 pub flags: [bool; 3],
130 pub origin: [f64; 3],
132}
133
134#[derive(Debug, Clone)]
136pub struct BoundarySpec {
137 pub elem_start: usize,
139 pub elem_end: usize,
141 pub bc_type: String,
143 pub value_re: f64,
145 pub curve_re: i32,
147 pub value_im: f64,
149 pub curve_im: i32,
151}
152
153#[derive(Debug, Clone)]
155pub struct PlaneWaveSource {
156 pub direction: [f64; 3],
158 pub amplitude_re: f64,
160 pub curve_re: i32,
162 pub amplitude_im: f64,
164 pub curve_im: i32,
166}
167
168#[derive(Debug, Clone)]
170pub struct PointSource {
171 pub position: [f64; 3],
173 pub amplitude_re: f64,
175 pub curve_re: i32,
177 pub amplitude_im: f64,
179 pub curve_im: i32,
181}
182
183#[derive(Debug, thiserror::Error)]
185pub enum NcParseError {
186 #[error("IO error: {0}")]
187 Io(#[from] std::io::Error),
188 #[error("Parse error at line {line}: {message}")]
189 Parse { line: usize, message: String },
190 #[error("Missing required section: {0}")]
191 MissingSection(String),
192 #[error("Invalid format: {0}")]
193 InvalidFormat(String),
194}
195
196pub fn parse_nc_input<P: AsRef<Path>>(path: P) -> Result<NcInputConfig, NcParseError> {
198 let path = path.as_ref();
199 let base_dir = path.parent().unwrap_or(Path::new(".")).to_path_buf();
200 let content = fs::read_to_string(path)?;
201
202 parse_nc_input_string(&content, base_dir)
203}
204
205pub fn parse_nc_input_string(content: &str, base_dir: PathBuf) -> Result<NcInputConfig, NcParseError> {
207 let mut config = NcInputConfig {
208 version: String::new(),
209 description: String::new(),
210 control_params_i: Vec::new(),
211 control_params_ii: Vec::new(),
212 frequency_curve: Vec::new(),
213 main_params_i: MainParamsI::default(),
214 main_params_ii: MainParamsII::default(),
215 main_params_iii: Vec::new(),
216 main_params_iv: MainParamsIV::default(),
217 node_files: Vec::new(),
218 element_files: Vec::new(),
219 symmetry: None,
220 boundary_conditions: Vec::new(),
221 plane_waves: Vec::new(),
222 point_sources: Vec::new(),
223 base_dir,
224 };
225
226 let lines: Vec<&str> = content.lines().collect();
227 let mut i = 0;
228
229 while i < lines.len() {
230 let line = lines[i].trim();
231
232 if line.is_empty() || line.starts_with("##") {
234 i += 1;
235 continue;
236 }
237
238 if line.starts_with('#') && !line.starts_with("##") {
240 i += 1;
242 continue;
243 }
244
245 if line.starts_with("Mesh2HRTF") {
247 config.version = line.to_string();
248 i += 1;
249 continue;
250 }
251
252 if config.version.is_empty() {
254 i += 1;
255 continue;
256 }
257
258 if i > 0 && lines.get(i.saturating_sub(1)).map_or(false, |l| {
260 l.contains("Controlparameter I") && !l.contains("Controlparameter II")
261 }) {
262 config.control_params_i = parse_int_line(line);
263 i += 1;
264 continue;
265 }
266
267 if i > 0 && lines.get(i.saturating_sub(1)).map_or(false, |l| l.contains("Controlparameter II")) {
269 config.control_params_ii = parse_float_line(line);
270 i += 1;
271 continue;
272 }
273
274 if i > 0 && lines.get(i.saturating_sub(1)).map_or(false, |l| l.contains("Frequency Curve")) {
276 let header = parse_int_line(line);
277 let num_points = header.get(1).copied().unwrap_or(0) as usize;
278 i += 1;
279
280 for _ in 0..num_points {
281 if i < lines.len() {
282 let values = parse_float_line(lines[i]);
283 if values.len() >= 3 {
284 config.frequency_curve.push((values[0], values[1], values[2]));
285 }
286 i += 1;
287 }
288 }
289 continue;
290 }
291
292 if i > 0 && lines.get(i.saturating_sub(1)).map_or(false, |l| {
294 l.contains("Main Parameters I") && !l.contains("Main Parameters II")
295 && !l.contains("Main Parameters III") && !l.contains("Main Parameters IV")
296 }) {
297 let values = parse_int_line(line);
298 config.main_params_i = MainParamsI {
299 element_type: values.get(0).copied().unwrap_or(0),
300 num_nodes: values.get(1).copied().unwrap_or(0) as usize,
301 num_elements: values.get(2).copied().unwrap_or(0) as usize,
302 num_object_files: values.get(3).copied().unwrap_or(0),
303 num_eval_files: values.get(4).copied().unwrap_or(0),
304 bc_type: values.get(5).copied().unwrap_or(0),
305 solver_method: values.get(6).copied().unwrap_or(0),
306 fmm_method: values.get(7).copied().unwrap_or(0),
307 parallel: values.get(8).copied().unwrap_or(0),
308 };
309 i += 1;
310 continue;
311 }
312
313 if i > 0 && lines.get(i.saturating_sub(1)).map_or(false, |l| {
315 l.contains("Main Parameters II") && !l.contains("Main Parameters III")
316 && !l.contains("Main Parameters IV")
317 }) {
318 let values = parse_mixed_line(line);
319 config.main_params_ii = MainParamsII {
320 preconditioner: values.get(0).map(|v| *v as i32).unwrap_or(0),
321 iterative_solver: values.get(1).map(|v| *v as i32).unwrap_or(0),
322 reserved1: values.get(2).map(|v| *v as i32).unwrap_or(0),
323 reserved2: values.get(3).copied().unwrap_or(0.0),
324 output_level: values.get(4).map(|v| *v as i32).unwrap_or(0),
325 reserved3: values.get(5).map(|v| *v as i32).unwrap_or(0),
326 reserved4: values.get(6).map(|v| *v as i32).unwrap_or(0),
327 };
328 i += 1;
329 continue;
330 }
331
332 if i > 0 && lines.get(i.saturating_sub(1)).map_or(false, |l| {
334 l.contains("Main Parameters III") && !l.contains("Main Parameters IV")
335 }) {
336 config.main_params_iii = parse_int_line(line);
337 i += 1;
338 continue;
339 }
340
341 if i > 0 && lines.get(i.saturating_sub(1)).map_or(false, |l| l.contains("Main Parameters IV")) {
343 let values = parse_float_line(line);
344 config.main_params_iv = MainParamsIV {
345 speed_of_sound: values.get(0).copied().unwrap_or(343.0),
346 density: values.get(1).copied().unwrap_or(1.21),
347 reference_pressure: values.get(2).copied().unwrap_or(1.0),
348 reserved: values.get(3..).map(|s| s.to_vec()).unwrap_or_default(),
349 };
350 i += 1;
351 continue;
352 }
353
354 if line == "NODES" {
356 i += 1;
357 while i < lines.len() {
358 let node_line = lines[i].trim();
359 if node_line.starts_with("##") || node_line.is_empty() {
360 break;
361 }
362 if !node_line.starts_with('#') {
363 let path = config.base_dir.join(node_line);
364 config.node_files.push(path);
365 }
366 i += 1;
367 }
368 continue;
369 }
370
371 if line == "ELEMENTS" {
373 i += 1;
374 while i < lines.len() {
375 let elem_line = lines[i].trim();
376 if elem_line.starts_with("##") || elem_line.is_empty() {
377 break;
378 }
379 if !elem_line.starts_with('#') {
380 let path = config.base_dir.join(elem_line);
381 config.element_files.push(path);
382 }
383 i += 1;
384 }
385 continue;
386 }
387
388 if line == "SYMMETRY" {
390 i += 1;
391 if i < lines.len() {
392 let flags_line = lines[i].trim();
393 if !flags_line.starts_with('#') {
394 let flags = parse_int_line(flags_line);
395 i += 1;
396 if i < lines.len() {
397 let origin = parse_float_line(lines[i].trim());
398 config.symmetry = Some(SymmetryConfig {
399 flags: [
400 flags.get(0).copied().unwrap_or(0) != 0,
401 flags.get(1).copied().unwrap_or(0) != 0,
402 flags.get(2).copied().unwrap_or(0) != 0,
403 ],
404 origin: [
405 origin.get(0).copied().unwrap_or(0.0),
406 origin.get(1).copied().unwrap_or(0.0),
407 origin.get(2).copied().unwrap_or(0.0),
408 ],
409 });
410 i += 1;
411 }
412 }
413 }
414 continue;
415 }
416
417 if line == "BOUNDARY" {
419 i += 1;
420 while i < lines.len() {
421 let bc_line = lines[i].trim();
422 if bc_line.starts_with("##") || bc_line == "RETU" {
423 i += 1;
424 break;
425 }
426 if bc_line.starts_with('#') || bc_line.is_empty() {
427 i += 1;
428 continue;
429 }
430
431 if let Some(bc) = parse_boundary_line(bc_line) {
432 config.boundary_conditions.push(bc);
433 }
434 i += 1;
435 }
436 continue;
437 }
438
439 if line == "PLANE WAVES" {
441 i += 1;
442 while i < lines.len() {
443 let pw_line = lines[i].trim();
444 if pw_line.starts_with("##") || pw_line.is_empty() {
445 break;
446 }
447 if !pw_line.starts_with('#') {
448 if let Some(pw) = parse_plane_wave_line(pw_line) {
449 config.plane_waves.push(pw);
450 }
451 }
452 i += 1;
453 }
454 continue;
455 }
456
457 if line == "POINT SOURCES" {
459 i += 1;
460 while i < lines.len() {
461 let ps_line = lines[i].trim();
462 if ps_line.starts_with("##") || ps_line.is_empty() {
463 break;
464 }
465 if !ps_line.starts_with('#') {
466 if let Some(ps) = parse_point_source_line(ps_line) {
467 config.point_sources.push(ps);
468 }
469 }
470 i += 1;
471 }
472 continue;
473 }
474
475 if line == "END" {
477 break;
478 }
479
480 i += 1;
481 }
482
483 Ok(config)
484}
485
486fn parse_int_line(line: &str) -> Vec<i32> {
488 line.split_whitespace()
489 .filter_map(|s| s.parse::<i32>().ok())
490 .collect()
491}
492
493fn parse_float_line(line: &str) -> Vec<f64> {
495 line.split_whitespace()
496 .filter_map(|s| {
497 let s = s.replace(" ", "");
499 s.parse::<f64>().ok()
500 })
501 .collect()
502}
503
504fn parse_mixed_line(line: &str) -> Vec<f64> {
506 parse_float_line(line)
507}
508
509fn parse_boundary_line(line: &str) -> Option<BoundarySpec> {
511 let parts: Vec<&str> = line.split_whitespace().collect();
512
513 if parts.len() >= 9 && parts[0] == "ELEM" && parts[2] == "TO" {
515 let elem_start = parts[1].parse().ok()?;
516 let elem_end = parts[3].parse().ok()?;
517 let bc_type = parts[4].to_string();
518 let value_re = parts[5].parse().ok()?;
519 let curve_re = parts[6].parse().ok()?;
520 let value_im = parts[7].parse().ok()?;
521 let curve_im = parts[8].parse().ok()?;
522
523 return Some(BoundarySpec {
524 elem_start,
525 elem_end,
526 bc_type,
527 value_re,
528 curve_re,
529 value_im,
530 curve_im,
531 });
532 }
533
534 None
535}
536
537fn parse_plane_wave_line(line: &str) -> Option<PlaneWaveSource> {
539 let values = parse_float_line(line);
540 if values.len() >= 8 {
541 Some(PlaneWaveSource {
542 direction: [values[1], values[2], values[3]],
543 amplitude_re: values[4],
544 curve_re: values[5] as i32,
545 amplitude_im: values[6],
546 curve_im: values[7] as i32,
547 })
548 } else {
549 None
550 }
551}
552
553fn parse_point_source_line(line: &str) -> Option<PointSource> {
555 let values = parse_float_line(line);
556 if values.len() >= 8 {
557 Some(PointSource {
558 position: [values[1], values[2], values[3]],
559 amplitude_re: values[4],
560 curve_re: values[5] as i32,
561 amplitude_im: values[6],
562 curve_im: values[7] as i32,
563 })
564 } else {
565 None
566 }
567}
568
569pub fn load_nc_nodes<P: AsRef<Path>>(path: P) -> Result<Array2<f64>, NcParseError> {
571 let content = fs::read_to_string(path)?;
572 let lines: Vec<&str> = content.lines().collect();
573
574 if lines.is_empty() {
575 return Ok(Array2::zeros((0, 3)));
576 }
577
578 let start_line = if lines[0].trim().parse::<usize>().is_ok() { 1 } else { 0 };
580
581 let mut nodes = Vec::new();
582 for line in &lines[start_line..] {
583 let values = parse_float_line(line);
584 if values.len() >= 4 {
585 nodes.push([values[1], values[2], values[3]]);
587 } else if values.len() >= 3 {
588 nodes.push([values[0], values[1], values[2]]);
590 }
591 }
592
593 let n = nodes.len();
594 let flat: Vec<f64> = nodes.into_iter().flatten().collect();
595 Array2::from_shape_vec((n, 3), flat)
596 .map_err(|e| NcParseError::InvalidFormat(e.to_string()))
597}
598
599pub fn load_nc_elements<P: AsRef<Path>>(
601 path: P,
602 property: ElementProperty,
603) -> Result<Vec<Element>, NcParseError> {
604 let content = fs::read_to_string(path)?;
605 let lines: Vec<&str> = content.lines().collect();
606
607 if lines.is_empty() {
608 return Ok(Vec::new());
609 }
610
611 let start_line = if lines[0].trim().parse::<usize>().is_ok() { 1 } else { 0 };
612
613 let mut elements = Vec::new();
614 for (idx, line) in lines[start_line..].iter().enumerate() {
615 let values: Vec<i32> = line
616 .split_whitespace()
617 .filter_map(|s| s.parse().ok())
618 .collect();
619
620 if values.len() >= 4 {
621 let connectivity: Vec<usize> = values[1..]
623 .iter()
624 .take_while(|&&v| v >= 0)
625 .map(|&v| v as usize)
626 .collect();
627
628 let element_type = if connectivity.len() == 3 {
629 ElementType::Tri3
630 } else {
631 ElementType::Quad4
632 };
633
634 let elem = Element {
635 connectivity,
636 element_type,
637 property,
638 normal: Array1::zeros(3),
639 node_normals: Array2::zeros((element_type.num_nodes(), 3)),
640 center: Array1::zeros(3),
641 area: 0.0,
642 boundary_condition: BoundaryCondition::Velocity(vec![Complex64::new(0.0, 0.0)]),
643 group: 0,
644 dof_addresses: vec![idx],
645 };
646 elements.push(elem);
647 }
648 }
649
650 Ok(elements)
651}
652
653impl NcInputConfig {
655 pub fn to_physics_params(&self, frequency: f64) -> PhysicsParams {
657 PhysicsParams::new(
658 frequency,
659 self.main_params_iv.speed_of_sound,
660 self.main_params_iv.density,
661 false, )
663 }
664
665 pub fn bem_method(&self) -> crate::core::types::BemMethod {
667 match self.main_params_i.solver_method {
668 0 => crate::core::types::BemMethod::Traditional,
669 1 => crate::core::types::BemMethod::SingleLevelFmm,
670 2 | 3 => crate::core::types::BemMethod::MultiLevelFmm,
671 _ => crate::core::types::BemMethod::Traditional,
672 }
673 }
674
675 pub fn solver_method(&self) -> crate::core::types::SolverMethod {
677 match self.main_params_ii.iterative_solver {
678 0 => crate::core::types::SolverMethod::Cgs,
679 1 => crate::core::types::SolverMethod::BiCgstab,
680 _ => crate::core::types::SolverMethod::Cgs,
681 }
682 }
683}
684
685#[cfg(test)]
686mod tests {
687 use super::*;
688
689 const SAMPLE_NC_INP: &str = r#"##-------------------------------------------
690## This file was created by mesh2input
691##-------------------------------------------
692Mesh2HRTF 1.0.0
693##
694Test Description
695##
696## Controlparameter I
6970 0 0 0 7 0
698##
699## Controlparameter II
7001 1 0.000001 0.00e+00 1 0 0
701##
702## Load Frequency Curve
7030 2
7040.000000 0.000000e+00 0.0
7050.000001 0.400000e+04 0.0
706##
707## 1. Main Parameters I
7082 100 50 0 0 2 1 0 0
709##
710## 2. Main Parameters II
7111 0 0 0.0000e+00 0 0 0
712##
713## 3. Main Parameters III
7140 0 0 0
715##
716## 4. Main Parameters IV
717343 1.21 1.0 0.0 0.0 0.0 0.0
718##
719NODES
720nodes.txt
721##
722ELEMENTS
723elements.txt
724##
725BOUNDARY
726ELEM 0 TO 49 VELO 1.0 -1 0.0 -1
727RETU
728##
729PLANE WAVES
7301 0.0 -1.0 0.0 1.0 -1 0.0 -1
731##
732END
733"#;
734
735 #[test]
736 fn test_parse_nc_input() {
737 let config = parse_nc_input_string(SAMPLE_NC_INP, PathBuf::from(".")).unwrap();
738
739 assert!(config.version.contains("Mesh2HRTF"));
740 assert_eq!(config.main_params_i.num_nodes, 100);
741 assert_eq!(config.main_params_i.num_elements, 50);
742 assert_eq!(config.main_params_i.solver_method, 1);
743 assert!((config.main_params_iv.speed_of_sound - 343.0).abs() < 0.01);
744 assert!((config.main_params_iv.density - 1.21).abs() < 0.01);
745 assert_eq!(config.node_files.len(), 1);
746 assert_eq!(config.element_files.len(), 1);
747 assert_eq!(config.boundary_conditions.len(), 1);
748 assert_eq!(config.plane_waves.len(), 1);
749 }
750
751 #[test]
752 fn test_parse_boundary_line() {
753 let bc = parse_boundary_line("ELEM 0 TO 100 VELO 1.0 -1 0.0 -1").unwrap();
754 assert_eq!(bc.elem_start, 0);
755 assert_eq!(bc.elem_end, 100);
756 assert_eq!(bc.bc_type, "VELO");
757 assert!((bc.value_re - 1.0).abs() < 0.001);
758 }
759
760 #[test]
761 fn test_parse_plane_wave() {
762 let pw = parse_plane_wave_line("1 0.0 -1.0 0.0 1.0 -1 0.0 -1").unwrap();
763 assert!((pw.direction[1] - (-1.0)).abs() < 0.001);
764 assert!((pw.amplitude_re - 1.0).abs() < 0.001);
765 }
766
767 #[test]
768 fn test_config_to_physics() {
769 let config = parse_nc_input_string(SAMPLE_NC_INP, PathBuf::from(".")).unwrap();
770 let physics = config.to_physics_params(1000.0);
771
772 assert!((physics.speed_of_sound - 343.0).abs() < 0.01);
773 assert!((physics.density - 1.21).abs() < 0.01);
774 assert!((physics.frequency - 1000.0).abs() < 0.01);
775 }
776}