1use crate::source::Source;
4use crate::types::{Point3D, RoomMesh, SurfaceElement};
5use serde::{Deserialize, Serialize};
6
7#[derive(Debug, Clone, Serialize, Deserialize)]
9pub enum RoomGeometry {
10 Rectangular(RectangularRoom),
12 LShaped(LShapedRoom),
14}
15
16impl RoomGeometry {
17 pub fn generate_mesh(&self, elements_per_meter: usize) -> RoomMesh {
19 match self {
20 RoomGeometry::Rectangular(room) => room.generate_mesh(elements_per_meter),
21 RoomGeometry::LShaped(room) => room.generate_mesh(elements_per_meter),
22 }
23 }
24
25 pub fn generate_adaptive_mesh(
27 &self,
28 base_elements_per_meter: usize,
29 frequency: f64,
30 sources: &[Source],
31 speed_of_sound: f64,
32 ) -> RoomMesh {
33 match self {
34 RoomGeometry::Rectangular(room) => room.generate_adaptive_mesh(
35 base_elements_per_meter,
36 frequency,
37 sources,
38 speed_of_sound,
39 ),
40 RoomGeometry::LShaped(room) => room.generate_adaptive_mesh(
41 base_elements_per_meter,
42 frequency,
43 sources,
44 speed_of_sound,
45 ),
46 }
47 }
48
49 pub fn get_edges(&self) -> Vec<(Point3D, Point3D)> {
51 match self {
52 RoomGeometry::Rectangular(room) => room.get_edges(),
53 RoomGeometry::LShaped(room) => room.get_edges(),
54 }
55 }
56
57 pub fn dimensions(&self) -> (f64, f64, f64) {
59 match self {
60 RoomGeometry::Rectangular(r) => (r.width, r.depth, r.height),
61 RoomGeometry::LShaped(r) => (r.width1.max(r.width2), r.depth1 + r.depth2, r.height),
62 }
63 }
64
65 pub fn volume(&self) -> f64 {
67 match self {
68 RoomGeometry::Rectangular(r) => r.width * r.depth * r.height,
69 RoomGeometry::LShaped(r) => {
70 r.width1 * r.depth1 * r.height + r.width2 * r.depth2 * r.height
71 }
72 }
73 }
74}
75
76struct SurfaceInfo {
78 origin: Point3D,
79 u_dir: Point3D,
80 v_dir: Point3D,
81 u_length: f64,
82 v_length: f64,
83}
84
85#[derive(Debug, Clone, Serialize, Deserialize)]
87pub struct RectangularRoom {
88 pub width: f64,
90 pub depth: f64,
92 pub height: f64,
94}
95
96impl RectangularRoom {
97 pub fn new(width: f64, depth: f64, height: f64) -> Self {
99 Self {
100 width,
101 depth,
102 height,
103 }
104 }
105
106 pub fn generate_mesh(&self, elements_per_meter: usize) -> RoomMesh {
108 let nx = (self.width * elements_per_meter as f64).ceil() as usize;
109 let ny = (self.depth * elements_per_meter as f64).ceil() as usize;
110 let nz = (self.height * elements_per_meter as f64).ceil() as usize;
111
112 let mut nodes = Vec::new();
113 let mut elements = Vec::new();
114
115 Self::add_surface_mesh(
117 &mut nodes,
118 &mut elements,
119 Point3D::new(0.0, 0.0, 0.0),
120 Point3D::new(self.width, 0.0, 0.0),
121 Point3D::new(0.0, self.depth, 0.0),
122 nx,
123 ny,
124 );
125
126 Self::add_surface_mesh(
128 &mut nodes,
129 &mut elements,
130 Point3D::new(0.0, 0.0, self.height),
131 Point3D::new(self.width, 0.0, self.height),
132 Point3D::new(0.0, self.depth, self.height),
133 nx,
134 ny,
135 );
136
137 Self::add_surface_mesh(
139 &mut nodes,
140 &mut elements,
141 Point3D::new(0.0, 0.0, 0.0),
142 Point3D::new(self.width, 0.0, 0.0),
143 Point3D::new(0.0, 0.0, self.height),
144 nx,
145 nz,
146 );
147
148 Self::add_surface_mesh(
150 &mut nodes,
151 &mut elements,
152 Point3D::new(0.0, self.depth, 0.0),
153 Point3D::new(self.width, self.depth, 0.0),
154 Point3D::new(0.0, self.depth, self.height),
155 nx,
156 nz,
157 );
158
159 Self::add_surface_mesh(
161 &mut nodes,
162 &mut elements,
163 Point3D::new(0.0, 0.0, 0.0),
164 Point3D::new(0.0, self.depth, 0.0),
165 Point3D::new(0.0, 0.0, self.height),
166 ny,
167 nz,
168 );
169
170 Self::add_surface_mesh(
172 &mut nodes,
173 &mut elements,
174 Point3D::new(self.width, 0.0, 0.0),
175 Point3D::new(self.width, self.depth, 0.0),
176 Point3D::new(self.width, 0.0, self.height),
177 ny,
178 nz,
179 );
180
181 RoomMesh { nodes, elements }
182 }
183
184 pub fn generate_adaptive_mesh(
186 &self,
187 base_elements_per_meter: usize,
188 frequency: f64,
189 sources: &[Source],
190 speed_of_sound: f64,
191 ) -> RoomMesh {
192 let wavelength = speed_of_sound / frequency;
193 let target_element_size = wavelength / 8.0;
194
195 let mut nodes = Vec::new();
196 let mut elements = Vec::new();
197
198 let surfaces = vec![
199 SurfaceInfo {
201 origin: Point3D::new(0.0, 0.0, 0.0),
202 u_dir: Point3D::new(self.width, 0.0, 0.0),
203 v_dir: Point3D::new(0.0, self.depth, 0.0),
204 u_length: self.width,
205 v_length: self.depth,
206 },
207 SurfaceInfo {
209 origin: Point3D::new(0.0, 0.0, self.height),
210 u_dir: Point3D::new(self.width, 0.0, self.height),
211 v_dir: Point3D::new(0.0, self.depth, self.height),
212 u_length: self.width,
213 v_length: self.depth,
214 },
215 SurfaceInfo {
217 origin: Point3D::new(0.0, 0.0, 0.0),
218 u_dir: Point3D::new(self.width, 0.0, 0.0),
219 v_dir: Point3D::new(0.0, 0.0, self.height),
220 u_length: self.width,
221 v_length: self.height,
222 },
223 SurfaceInfo {
225 origin: Point3D::new(0.0, self.depth, 0.0),
226 u_dir: Point3D::new(self.width, self.depth, 0.0),
227 v_dir: Point3D::new(0.0, self.depth, self.height),
228 u_length: self.width,
229 v_length: self.height,
230 },
231 SurfaceInfo {
233 origin: Point3D::new(0.0, 0.0, 0.0),
234 u_dir: Point3D::new(0.0, self.depth, 0.0),
235 v_dir: Point3D::new(0.0, 0.0, self.height),
236 u_length: self.depth,
237 v_length: self.height,
238 },
239 SurfaceInfo {
241 origin: Point3D::new(self.width, 0.0, 0.0),
242 u_dir: Point3D::new(self.width, self.depth, 0.0),
243 v_dir: Point3D::new(self.width, 0.0, self.height),
244 u_length: self.depth,
245 v_length: self.height,
246 },
247 ];
248
249 for surface in surfaces {
250 self.add_adaptive_surface_mesh(
251 &mut nodes,
252 &mut elements,
253 &surface,
254 target_element_size,
255 base_elements_per_meter,
256 sources,
257 );
258 }
259
260 RoomMesh { nodes, elements }
261 }
262
263 fn add_adaptive_surface_mesh(
264 &self,
265 nodes: &mut Vec<Point3D>,
266 elements: &mut Vec<SurfaceElement>,
267 surface: &SurfaceInfo,
268 target_element_size: f64,
269 base_elements_per_meter: usize,
270 sources: &[Source],
271 ) {
272 let nu_base = (surface.u_length / target_element_size).ceil() as usize;
273 let nv_base = (surface.v_length / target_element_size).ceil() as usize;
274
275 let nu_min = (surface.u_length * base_elements_per_meter as f64).ceil() as usize;
276 let nv_min = (surface.v_length * base_elements_per_meter as f64).ceil() as usize;
277
278 let nu = nu_base.max(nu_min);
279 let nv = nv_base.max(nv_min);
280
281 let near_source = sources.iter().any(|source| {
282 let dist_to_surface = self.distance_point_to_surface(&source.position, surface);
283 dist_to_surface < target_element_size * 2.0
284 });
285
286 let (nu_final, nv_final) = if near_source {
287 (nu * 2, nv * 2)
288 } else {
289 (nu, nv)
290 };
291
292 let base_idx = nodes.len();
293
294 for j in 0..=nv_final {
295 for i in 0..=nu_final {
296 let u = Self::graded_parameter(i as f64 / nu_final as f64);
297 let v = Self::graded_parameter(j as f64 / nv_final as f64);
298
299 let u_vec = Point3D::new(
300 surface.u_dir.x - surface.origin.x,
301 surface.u_dir.y - surface.origin.y,
302 surface.u_dir.z - surface.origin.z,
303 );
304 let v_vec = Point3D::new(
305 surface.v_dir.x - surface.origin.x,
306 surface.v_dir.y - surface.origin.y,
307 surface.v_dir.z - surface.origin.z,
308 );
309
310 let node = Point3D::new(
311 surface.origin.x + u * u_vec.x + v * v_vec.x,
312 surface.origin.y + u * u_vec.y + v * v_vec.y,
313 surface.origin.z + u * u_vec.z + v * v_vec.z,
314 );
315
316 nodes.push(node);
317 }
318 }
319
320 for j in 0..nv_final {
321 for i in 0..nu_final {
322 let i0 = base_idx + j * (nu_final + 1) + i;
323 let i1 = i0 + 1;
324 let i2 = i0 + (nu_final + 1);
325 let i3 = i2 + 1;
326
327 elements.push(SurfaceElement::triangle(i0, i1, i2));
328 elements.push(SurfaceElement::triangle(i1, i3, i2));
329 }
330 }
331 }
332
333 fn graded_parameter(t: f64) -> f64 {
334 let grading = 0.15;
335
336 if t < grading {
337 0.5 * (t / grading).powi(2) * grading
338 } else if t > 1.0 - grading {
339 let t_rel = (t - (1.0 - grading)) / grading;
340 1.0 - 0.5 * grading * (1.0 - t_rel).powi(2)
341 } else {
342 let t_mid = (t - grading) / (1.0 - 2.0 * grading);
343 grading + t_mid * (1.0 - 2.0 * grading)
344 }
345 }
346
347 fn distance_point_to_surface(&self, point: &Point3D, surface: &SurfaceInfo) -> f64 {
348 let u_vec = Point3D::new(
349 surface.u_dir.x - surface.origin.x,
350 surface.u_dir.y - surface.origin.y,
351 surface.u_dir.z - surface.origin.z,
352 );
353 let v_vec = Point3D::new(
354 surface.v_dir.x - surface.origin.x,
355 surface.v_dir.y - surface.origin.y,
356 surface.v_dir.z - surface.origin.z,
357 );
358
359 let normal = u_vec.cross(&v_vec);
360 let normal_length = normal.length();
361
362 if normal_length < 1e-10 {
363 return point.distance_to(&surface.origin);
364 }
365
366 let nx = normal.x / normal_length;
367 let ny = normal.y / normal_length;
368 let nz = normal.z / normal_length;
369
370 let dx = point.x - surface.origin.x;
371 let dy = point.y - surface.origin.y;
372 let dz = point.z - surface.origin.z;
373
374 (dx * nx + dy * ny + dz * nz).abs()
375 }
376
377 pub fn get_edges(&self) -> Vec<(Point3D, Point3D)> {
379 vec![
380 (
382 Point3D::new(0.0, 0.0, 0.0),
383 Point3D::new(self.width, 0.0, 0.0),
384 ),
385 (
386 Point3D::new(self.width, 0.0, 0.0),
387 Point3D::new(self.width, self.depth, 0.0),
388 ),
389 (
390 Point3D::new(self.width, self.depth, 0.0),
391 Point3D::new(0.0, self.depth, 0.0),
392 ),
393 (
394 Point3D::new(0.0, self.depth, 0.0),
395 Point3D::new(0.0, 0.0, 0.0),
396 ),
397 (
399 Point3D::new(0.0, 0.0, self.height),
400 Point3D::new(self.width, 0.0, self.height),
401 ),
402 (
403 Point3D::new(self.width, 0.0, self.height),
404 Point3D::new(self.width, self.depth, self.height),
405 ),
406 (
407 Point3D::new(self.width, self.depth, self.height),
408 Point3D::new(0.0, self.depth, self.height),
409 ),
410 (
411 Point3D::new(0.0, self.depth, self.height),
412 Point3D::new(0.0, 0.0, self.height),
413 ),
414 (
416 Point3D::new(0.0, 0.0, 0.0),
417 Point3D::new(0.0, 0.0, self.height),
418 ),
419 (
420 Point3D::new(self.width, 0.0, 0.0),
421 Point3D::new(self.width, 0.0, self.height),
422 ),
423 (
424 Point3D::new(self.width, self.depth, 0.0),
425 Point3D::new(self.width, self.depth, self.height),
426 ),
427 (
428 Point3D::new(0.0, self.depth, 0.0),
429 Point3D::new(0.0, self.depth, self.height),
430 ),
431 ]
432 }
433
434 fn add_surface_mesh(
435 nodes: &mut Vec<Point3D>,
436 elements: &mut Vec<SurfaceElement>,
437 origin: Point3D,
438 u_dir: Point3D,
439 v_dir: Point3D,
440 nu: usize,
441 nv: usize,
442 ) {
443 let base_idx = nodes.len();
444
445 for j in 0..=nv {
446 for i in 0..=nu {
447 let u = i as f64 / nu as f64;
448 let v = j as f64 / nv as f64;
449
450 let node = Point3D::new(
451 origin.x + u * (u_dir.x - origin.x) + v * (v_dir.x - origin.x),
452 origin.y + u * (u_dir.y - origin.y) + v * (v_dir.y - origin.y),
453 origin.z + u * (u_dir.z - origin.z) + v * (v_dir.z - origin.z),
454 );
455 nodes.push(node);
456 }
457 }
458
459 for j in 0..nv {
460 for i in 0..nu {
461 let n0 = base_idx + j * (nu + 1) + i;
462 let n1 = base_idx + j * (nu + 1) + i + 1;
463 let n2 = base_idx + (j + 1) * (nu + 1) + i + 1;
464 let n3 = base_idx + (j + 1) * (nu + 1) + i;
465
466 elements.push(SurfaceElement::quad(n0, n1, n2, n3));
467 }
468 }
469 }
470}
471
472#[derive(Debug, Clone, Serialize, Deserialize)]
474pub struct LShapedRoom {
475 pub width1: f64,
477 pub depth1: f64,
479 pub width2: f64,
481 pub depth2: f64,
483 pub height: f64,
485}
486
487impl LShapedRoom {
488 pub fn new(width1: f64, depth1: f64, width2: f64, depth2: f64, height: f64) -> Self {
490 Self {
491 width1,
492 depth1,
493 width2,
494 depth2,
495 height,
496 }
497 }
498
499 pub fn generate_mesh(&self, elements_per_meter: usize) -> RoomMesh {
501 let mut nodes = Vec::new();
502 let mut elements = Vec::new();
503
504 let nx1 = (self.width1 * elements_per_meter as f64).ceil() as usize;
505 let ny1 = (self.depth1 * elements_per_meter as f64).ceil() as usize;
506 let nx2 = (self.width2 * elements_per_meter as f64).ceil() as usize;
507 let ny2 = (self.depth2 * elements_per_meter as f64).ceil() as usize;
508 let nz = (self.height * elements_per_meter as f64).ceil() as usize;
509
510 Self::add_surface_mesh_lshaped(
512 &mut nodes,
513 &mut elements,
514 Point3D::new(0.0, 0.0, 0.0),
515 Point3D::new(self.width1, 0.0, 0.0),
516 Point3D::new(0.0, self.depth1, 0.0),
517 nx1,
518 ny1,
519 );
520
521 Self::add_surface_mesh_lshaped(
523 &mut nodes,
524 &mut elements,
525 Point3D::new(0.0, self.depth1, 0.0),
526 Point3D::new(self.width2, self.depth1, 0.0),
527 Point3D::new(0.0, self.depth1 + self.depth2, 0.0),
528 nx2,
529 ny2,
530 );
531
532 Self::add_surface_mesh_lshaped(
534 &mut nodes,
535 &mut elements,
536 Point3D::new(0.0, 0.0, self.height),
537 Point3D::new(self.width1, 0.0, self.height),
538 Point3D::new(0.0, self.depth1, self.height),
539 nx1,
540 ny1,
541 );
542
543 Self::add_surface_mesh_lshaped(
545 &mut nodes,
546 &mut elements,
547 Point3D::new(0.0, self.depth1, self.height),
548 Point3D::new(self.width2, self.depth1, self.height),
549 Point3D::new(0.0, self.depth1 + self.depth2, self.height),
550 nx2,
551 ny2,
552 );
553
554 Self::add_surface_mesh_lshaped(
557 &mut nodes,
558 &mut elements,
559 Point3D::new(0.0, 0.0, 0.0),
560 Point3D::new(self.width1, 0.0, 0.0),
561 Point3D::new(0.0, 0.0, self.height),
562 nx1,
563 nz,
564 );
565
566 Self::add_surface_mesh_lshaped(
568 &mut nodes,
569 &mut elements,
570 Point3D::new(self.width1, 0.0, 0.0),
571 Point3D::new(self.width1, self.depth1, 0.0),
572 Point3D::new(self.width1, 0.0, self.height),
573 ny1,
574 nz,
575 );
576
577 let total_depth = self.depth1 + self.depth2;
579 let ny_total = ((total_depth) * elements_per_meter as f64).ceil() as usize;
580 Self::add_surface_mesh_lshaped(
581 &mut nodes,
582 &mut elements,
583 Point3D::new(0.0, 0.0, 0.0),
584 Point3D::new(0.0, total_depth, 0.0),
585 Point3D::new(0.0, 0.0, self.height),
586 ny_total,
587 nz,
588 );
589
590 Self::add_surface_mesh_lshaped(
592 &mut nodes,
593 &mut elements,
594 Point3D::new(0.0, total_depth, 0.0),
595 Point3D::new(self.width2, total_depth, 0.0),
596 Point3D::new(0.0, total_depth, self.height),
597 nx2,
598 nz,
599 );
600
601 Self::add_surface_mesh_lshaped(
603 &mut nodes,
604 &mut elements,
605 Point3D::new(self.width2, self.depth1, 0.0),
606 Point3D::new(self.width2, total_depth, 0.0),
607 Point3D::new(self.width2, self.depth1, self.height),
608 ny2,
609 nz,
610 );
611
612 let internal_width = self.width1 - self.width2;
614 let nx_internal = (internal_width * elements_per_meter as f64).ceil() as usize;
615 Self::add_surface_mesh_lshaped(
616 &mut nodes,
617 &mut elements,
618 Point3D::new(self.width2, self.depth1, 0.0),
619 Point3D::new(self.width1, self.depth1, 0.0),
620 Point3D::new(self.width2, self.depth1, self.height),
621 nx_internal,
622 nz,
623 );
624
625 RoomMesh { nodes, elements }
626 }
627
628 pub fn get_edges(&self) -> Vec<(Point3D, Point3D)> {
630 let total_depth = self.depth1 + self.depth2;
631 vec![
632 (
634 Point3D::new(0.0, 0.0, 0.0),
635 Point3D::new(self.width1, 0.0, 0.0),
636 ),
637 (
638 Point3D::new(self.width1, 0.0, 0.0),
639 Point3D::new(self.width1, self.depth1, 0.0),
640 ),
641 (
642 Point3D::new(self.width1, self.depth1, 0.0),
643 Point3D::new(self.width2, self.depth1, 0.0),
644 ),
645 (
646 Point3D::new(self.width2, self.depth1, 0.0),
647 Point3D::new(self.width2, total_depth, 0.0),
648 ),
649 (
650 Point3D::new(self.width2, total_depth, 0.0),
651 Point3D::new(0.0, total_depth, 0.0),
652 ),
653 (
654 Point3D::new(0.0, total_depth, 0.0),
655 Point3D::new(0.0, 0.0, 0.0),
656 ),
657 (
659 Point3D::new(0.0, 0.0, self.height),
660 Point3D::new(self.width1, 0.0, self.height),
661 ),
662 (
663 Point3D::new(self.width1, 0.0, self.height),
664 Point3D::new(self.width1, self.depth1, self.height),
665 ),
666 (
667 Point3D::new(self.width1, self.depth1, self.height),
668 Point3D::new(self.width2, self.depth1, self.height),
669 ),
670 (
671 Point3D::new(self.width2, self.depth1, self.height),
672 Point3D::new(self.width2, total_depth, self.height),
673 ),
674 (
675 Point3D::new(self.width2, total_depth, self.height),
676 Point3D::new(0.0, total_depth, self.height),
677 ),
678 (
679 Point3D::new(0.0, total_depth, self.height),
680 Point3D::new(0.0, 0.0, self.height),
681 ),
682 (
684 Point3D::new(0.0, 0.0, 0.0),
685 Point3D::new(0.0, 0.0, self.height),
686 ),
687 (
688 Point3D::new(self.width1, 0.0, 0.0),
689 Point3D::new(self.width1, 0.0, self.height),
690 ),
691 (
692 Point3D::new(self.width1, self.depth1, 0.0),
693 Point3D::new(self.width1, self.depth1, self.height),
694 ),
695 (
696 Point3D::new(self.width2, self.depth1, 0.0),
697 Point3D::new(self.width2, self.depth1, self.height),
698 ),
699 (
700 Point3D::new(self.width2, total_depth, 0.0),
701 Point3D::new(self.width2, total_depth, self.height),
702 ),
703 (
704 Point3D::new(0.0, total_depth, 0.0),
705 Point3D::new(0.0, total_depth, self.height),
706 ),
707 ]
708 }
709
710 fn add_surface_mesh_lshaped(
711 nodes: &mut Vec<Point3D>,
712 elements: &mut Vec<SurfaceElement>,
713 origin: Point3D,
714 u_dir: Point3D,
715 v_dir: Point3D,
716 nu: usize,
717 nv: usize,
718 ) {
719 let base_idx = nodes.len();
720
721 for j in 0..=nv {
722 for i in 0..=nu {
723 let u = i as f64 / nu as f64;
724 let v = j as f64 / nv as f64;
725
726 let node = Point3D::new(
727 origin.x + u * (u_dir.x - origin.x) + v * (v_dir.x - origin.x),
728 origin.y + u * (u_dir.y - origin.y) + v * (v_dir.y - origin.y),
729 origin.z + u * (u_dir.z - origin.z) + v * (v_dir.z - origin.z),
730 );
731 nodes.push(node);
732 }
733 }
734
735 for j in 0..nv {
736 for i in 0..nu {
737 let n0 = base_idx + j * (nu + 1) + i;
738 let n1 = base_idx + j * (nu + 1) + i + 1;
739 let n2 = base_idx + (j + 1) * (nu + 1) + i + 1;
740 let n3 = base_idx + (j + 1) * (nu + 1) + i;
741
742 elements.push(SurfaceElement::quad(n0, n1, n2, n3));
743 }
744 }
745 }
746
747 pub fn generate_adaptive_mesh(
749 &self,
750 base_elements_per_meter: usize,
751 _frequency: f64,
752 _sources: &[Source],
753 _speed_of_sound: f64,
754 ) -> RoomMesh {
755 self.generate_mesh(base_elements_per_meter)
757 }
758}
759
760#[cfg(test)]
761mod tests {
762 use super::*;
763
764 #[test]
765 fn test_rectangular_room() {
766 let room = RectangularRoom::new(5.0, 4.0, 3.0);
767 assert_eq!(room.width, 5.0);
768 assert_eq!(room.depth, 4.0);
769 assert_eq!(room.height, 3.0);
770 }
771
772 #[test]
773 fn test_rectangular_room_mesh() {
774 let room = RectangularRoom::new(2.0, 2.0, 2.0);
775 let mesh = room.generate_mesh(1);
776 assert!(mesh.num_nodes() > 0);
777 assert!(mesh.num_elements() > 0);
778 }
779
780 #[test]
781 fn test_lshaped_room() {
782 let room = LShapedRoom::new(5.0, 4.0, 3.0, 3.0, 2.5);
783 assert_eq!(room.width1, 5.0);
784 assert_eq!(room.depth1, 4.0);
785 assert_eq!(room.width2, 3.0);
786 assert_eq!(room.depth2, 3.0);
787 assert_eq!(room.height, 2.5);
788 }
789}