1use crate::{
4 grid::local_grid::{SingleElementGrid, SingleElementGridBuilder},
5 traits::Builder,
6 types::Scalar,
7};
8use ndelement::{ciarlet::CiarletElement, map::IdentityMap, types::ReferenceCellType};
9
10pub fn unit_interval<T: Scalar>(
14 nx: usize,
15) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
16 let mut b = SingleElementGridBuilder::new_with_capacity(
17 1,
18 nx + 1,
19 nx,
20 (ReferenceCellType::Interval, 1),
21 );
22 for i in 0..nx + 1 {
23 b.add_point(i, &[T::from(i).unwrap() / T::from(nx).unwrap()]);
24 }
25
26 for i in 0..nx {
27 b.add_cell(i, &[i, i + 1]);
28 }
29
30 b.create_grid()
31}
32
33pub fn unit_square<T: Scalar>(
37 nx: usize,
38 ny: usize,
39 cell_type: ReferenceCellType,
40) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
41 let mut b = SingleElementGridBuilder::new_with_capacity(
42 2,
43 (nx + 1) * (ny + 1),
44 match cell_type {
45 ReferenceCellType::Triangle => 2 * nx * ny,
46 ReferenceCellType::Quadrilateral => 2 * nx * ny,
47 _ => {
48 panic!("Unsupported cell type: {cell_type:?}")
49 }
50 },
51 (cell_type, 1),
52 );
53 for i in 0..nx + 1 {
54 for j in 0..ny + 1 {
55 b.add_point(
56 i * (ny + 1) + j,
57 &[
58 T::from(i).unwrap() / T::from(nx).unwrap(),
59 T::from(j).unwrap() / T::from(ny).unwrap(),
60 ],
61 );
62 }
63 }
64
65 let dx = ny + 1;
66 let dy = 1;
67 match cell_type {
68 ReferenceCellType::Triangle => {
69 for i in 0..nx {
70 for j in 0..ny {
71 let origin = i * dx + j * dy;
72 b.add_cell(2 * (j * nx + i), &[origin, origin + dx, origin + dx + dy]);
73 b.add_cell(
74 2 * (j * nx + i) + 1,
75 &[origin, origin + dx + dy, origin + dy],
76 );
77 }
78 }
79 }
80 ReferenceCellType::Quadrilateral => {
81 for i in 0..nx {
82 for j in 0..ny {
83 let origin = i * dx + j * dy;
84 b.add_cell(
85 j * nx + i,
86 &[origin, origin + dx, origin + dy, origin + dx + dy],
87 );
88 }
89 }
90 }
91 _ => {
92 panic!("Unsupported cell type: {cell_type:?}")
93 }
94 }
95
96 b.create_grid()
97}
98
99pub fn unit_square_boundary<T: Scalar>(
103 nx: usize,
104 ny: usize,
105) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
106 let mut b = SingleElementGridBuilder::new_with_capacity(
107 2,
108 2 * (nx + ny),
109 2 * (nx + ny),
110 (ReferenceCellType::Interval, 1),
111 );
112 let dx = ny + 1;
113 let dy = 1;
114
115 for i in 0..nx + 1 {
116 b.add_point(
117 i * dx,
118 &[T::from(i).unwrap() / T::from(nx).unwrap(), T::zero()],
119 );
120 b.add_point(
121 i * dx + ny * dy,
122 &[T::from(i).unwrap() / T::from(nx).unwrap(), T::one()],
123 );
124 }
125 for j in 1..ny {
126 b.add_point(
127 j * dy,
128 &[T::zero(), T::from(j).unwrap() / T::from(ny).unwrap()],
129 );
130 b.add_point(
131 nx * dx + j * dy,
132 &[T::one(), T::from(j).unwrap() / T::from(ny).unwrap()],
133 );
134 }
135
136 for i in 0..nx {
137 let origin = i * dx;
138 b.add_cell(2 * i, &[origin, origin + dx]);
139 let origin = i * dx + ny * dy;
140 b.add_cell(2 * i + 1, &[origin, origin + dx]);
141 }
142 for j in 0..ny {
143 let origin = j * dy;
144 b.add_cell(2 * nx + 2 * j, &[origin, origin + dy]);
145 let origin = nx * dx + j * dy;
146 b.add_cell(2 * nx + 2 * j + 1, &[origin, origin + dy]);
147 }
148
149 b.create_grid()
150}
151
152pub fn unit_cube<T: Scalar>(
157 nx: usize,
158 ny: usize,
159 nz: usize,
160 cell_type: ReferenceCellType,
161) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
162 let mut b = SingleElementGridBuilder::new_with_capacity(
163 3,
164 (nx + 1) * (ny + 1) * (nz + 1),
165 match cell_type {
166 ReferenceCellType::Tetrahedron => 6 * nx * ny * nz,
167 ReferenceCellType::Hexahedron => 2 * nx * ny * nz,
168 _ => {
169 panic!("Unsupported cell type: {cell_type:?}")
170 }
171 },
172 (cell_type, 1),
173 );
174 for i in 0..nx + 1 {
175 for j in 0..ny + 1 {
176 for k in 0..nz + 1 {
177 b.add_point(
178 (i * (ny + 1) + j) * (nz + 1) + k,
179 &[
180 T::from(i).unwrap() / T::from(nx).unwrap(),
181 T::from(j).unwrap() / T::from(ny).unwrap(),
182 T::from(k).unwrap() / T::from(nz).unwrap(),
183 ],
184 );
185 }
186 }
187 }
188
189 let dx = (ny + 1) * (nz + 1);
190 let dy = nz + 1;
191 let dz = 1;
192 match cell_type {
193 ReferenceCellType::Tetrahedron => {
194 for i in 0..nx {
195 for j in 0..ny {
196 for k in 0..nz {
197 let origin = i * dx + j * dy + k * dz;
198 b.add_cell(
199 6 * ((j * nx + i) * ny + k),
200 &[origin, origin + dx, origin + dx + dy, origin + dx + dy + dz],
201 );
202 b.add_cell(
203 6 * ((j * nx + i) * ny + k) + 1,
204 &[origin, origin + dy, origin + dx + dy, origin + dx + dy + dz],
205 );
206 b.add_cell(
207 6 * ((j * nx + i) * ny + k) + 2,
208 &[origin, origin + dx, origin + dx + dz, origin + dx + dy + dz],
209 );
210 b.add_cell(
211 6 * ((j * nx + i) * ny + k) + 3,
212 &[origin, origin + dz, origin + dx + dz, origin + dx + dy + dz],
213 );
214 b.add_cell(
215 6 * ((j * nx + i) * ny + k) + 4,
216 &[origin, origin + dy, origin + dy + dz, origin + dx + dy + dz],
217 );
218 b.add_cell(
219 6 * ((j * nx + i) * ny + k) + 5,
220 &[origin, origin + dz, origin + dy + dz, origin + dx + dy + dz],
221 );
222 }
223 }
224 }
225 }
226 ReferenceCellType::Hexahedron => {
227 for i in 0..nx {
228 for j in 0..ny {
229 for k in 0..nz {
230 let origin = i * dx + j * dy + k * dz;
231 b.add_cell(
232 (j * nx + i) * ny + k,
233 &[
234 origin,
235 origin + dx,
236 origin + dy,
237 origin + dy,
238 origin + dz,
239 origin + dx + dz,
240 origin + dy + dz,
241 origin + dy + dz,
242 ],
243 );
244 }
245 }
246 }
247 }
248 _ => {
249 panic!("Unsupported cell type: {cell_type:?}")
250 }
251 }
252
253 b.create_grid()
254}
255
256pub fn unit_cube_boundary<T: Scalar>(
261 nx: usize,
262 ny: usize,
263 nz: usize,
264 cell_type: ReferenceCellType,
265) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
266 let mut b = SingleElementGridBuilder::new_with_capacity(
267 3,
268 (nx + 1) * (ny + 1) * (nz + 1) - (nx - 1) * (ny - 1) * (nz - 1),
269 match cell_type {
270 ReferenceCellType::Triangle => 4 * (nx * ny + nx * nz + ny * nz),
271 ReferenceCellType::Quadrilateral => 2 * (nx * ny + nx * nz + ny * nz),
272 _ => {
273 panic!("Unsupported cell type: {cell_type:?}")
274 }
275 },
276 (cell_type, 1),
277 );
278 for i in 0..nx + 1 {
279 for j in 0..ny + 1 {
280 for k in if i == 0 || i == nx || j == 0 || j == ny {
281 (0..nz + 1).collect::<Vec<_>>()
282 } else {
283 vec![0, nz]
284 } {
285 b.add_point(
286 (i * (ny + 1) + j) * (nz + 1) + k,
287 &[
288 T::from(i).unwrap() / T::from(nx).unwrap(),
289 T::from(j).unwrap() / T::from(ny).unwrap(),
290 T::from(k).unwrap() / T::from(nz).unwrap(),
291 ],
292 );
293 }
294 }
295 }
296
297 let dx = (ny + 1) * (nz + 1);
298 let dy = nz + 1;
299 let dz = 1;
300 match cell_type {
301 ReferenceCellType::Triangle => {
302 let mut cell_n = 0;
303 for i in 0..nx {
304 for j in 0..ny {
305 let origin = i * dx + j * dy;
306 b.add_cell(cell_n, &[origin, origin + dx, origin + dx + dy]);
307 b.add_cell(cell_n + 1, &[origin, origin + dx + dy, origin + dy]);
308 let origin = i * dx + j * dy + nz * dz;
309 b.add_cell(cell_n + 2, &[origin, origin + dx, origin + dx + dy]);
310 b.add_cell(cell_n + 3, &[origin, origin + dx + dy, origin + dy]);
311 cell_n += 4;
312 }
313 }
314 for i in 0..nx {
315 for k in 0..nz {
316 let origin = i * dx + k * dz;
317 b.add_cell(cell_n, &[origin, origin + dx, origin + dx + dz]);
318 b.add_cell(cell_n + 1, &[origin, origin + dx + dz, origin + dz]);
319 let origin = i * dx + ny * dy + k * dz;
320 b.add_cell(cell_n + 2, &[origin, origin + dx, origin + dx + dz]);
321 b.add_cell(cell_n + 3, &[origin, origin + dx + dz, origin + dz]);
322 cell_n += 4;
323 }
324 }
325 for j in 0..ny {
326 for k in 0..nz {
327 let origin = j * dy + k * dz;
328 b.add_cell(cell_n, &[origin, origin + dy, origin + dy + dz]);
329 b.add_cell(cell_n + 1, &[origin, origin + dy + dz, origin + dz]);
330 let origin = nx * dx + j * dy + k * dz;
331 b.add_cell(cell_n + 2, &[origin, origin + dy, origin + dy + dz]);
332 b.add_cell(cell_n + 3, &[origin, origin + dy + dz, origin + dz]);
333 cell_n += 4;
334 }
335 }
336 }
337 ReferenceCellType::Quadrilateral => {
338 let mut cell_n = 0;
339 for i in 0..nx {
340 for j in 0..ny {
341 let origin = i * dx + j * dy;
342 b.add_cell(
343 cell_n,
344 &[origin, origin + dx, origin + dy, origin + dx + dy],
345 );
346 let origin = i * dx + j * dy + nz * dz;
347 b.add_cell(
348 cell_n + 1,
349 &[origin, origin + dx, origin + dy, origin + dx + dy],
350 );
351 cell_n += 2;
352 }
353 }
354 for i in 0..nx {
355 for k in 0..nz {
356 let origin = i * dx + k * dz;
357 b.add_cell(
358 cell_n,
359 &[origin, origin + dx, origin + dz, origin + dx + dz],
360 );
361 let origin = i * dx + ny * dy + k * dz;
362 b.add_cell(
363 cell_n + 1,
364 &[origin, origin + dx, origin + dz, origin + dx + dz],
365 );
366 cell_n += 2;
367 }
368 }
369 for j in 0..ny {
370 for k in 0..nz {
371 let origin = j * dy + k * dz;
372 b.add_cell(
373 cell_n,
374 &[origin, origin + dy, origin + dz, origin + dy + dz],
375 );
376 let origin = nx * dx + j * dy + k * dz;
377 b.add_cell(
378 cell_n + 1,
379 &[origin, origin + dy, origin + dz, origin + dy + dz],
380 );
381 cell_n += 2;
382 }
383 }
384 }
385 _ => {
386 panic!("Unsupported cell type: {cell_type:?}")
387 }
388 }
389
390 b.create_grid()
391}
392
393pub fn unit_cube_edges<T: Scalar>(
398 nx: usize,
399 ny: usize,
400 nz: usize,
401) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
402 let mut b = SingleElementGridBuilder::new_with_capacity(
403 3,
404 4 * (nx + ny + nz + 1),
405 4 * (nx + ny + nz),
406 (ReferenceCellType::Interval, 1),
407 );
408 for i in 0..nx + 1 {
409 for j in if i == 0 || i == nx {
410 (0..ny + 1).collect::<Vec<_>>()
411 } else {
412 vec![0, ny]
413 } {
414 for k in if (i == 0 || i == nx) && (j == 0 || j == ny) {
415 (0..nz + 1).collect::<Vec<_>>()
416 } else {
417 vec![0, nz]
418 } {
419 b.add_point(
420 (i * (ny + 1) + j) * (nz + 1) + k,
421 &[
422 T::from(i).unwrap() / T::from(nx).unwrap(),
423 T::from(j).unwrap() / T::from(ny).unwrap(),
424 T::from(k).unwrap() / T::from(nz).unwrap(),
425 ],
426 );
427 }
428 }
429 }
430
431 let dx = (ny + 1) * (nz + 1);
432 let dy = nz + 1;
433 let dz = 1;
434 let mut cell_n = 0;
435 for i in 0..nx {
436 for origin in [
437 i * dx,
438 i * dx + ny * dy,
439 i * dx + nz * dz,
440 i * dx + ny * dy + nz * dz,
441 ] {
442 b.add_cell(cell_n, &[origin, origin + dx]);
443 cell_n += 1;
444 }
445 }
446 for j in 0..ny {
447 for origin in [
448 j * dy,
449 j * dy + nx * dx,
450 j * dy + nz * dz,
451 j * dy + nx * dx + nz * dz,
452 ] {
453 b.add_cell(cell_n, &[origin, origin + dy]);
454 cell_n += 1;
455 }
456 }
457 for k in 0..nz {
458 for origin in [
459 k * dz,
460 k * dz + nx * dx,
461 k * dz + ny * dy,
462 k * dz + nx * dx + ny * dy,
463 ] {
464 b.add_cell(cell_n, &[origin, origin + dz]);
465 cell_n += 1;
466 }
467 }
468
469 b.create_grid()
470}
471
472#[cfg(test)]
473mod test {
474 use super::*;
475 use crate::traits::{Entity, Geometry, Grid, Point};
476 use approx::*;
477 use itertools::izip;
478
479 fn max(values: &[f64]) -> f64 {
480 let mut out = values[0];
481 for i in &values[1..] {
482 if *i > out {
483 out = *i;
484 }
485 }
486 out
487 }
488
489 fn check_volume(
490 grid: &impl Grid<T = f64, EntityDescriptor = ReferenceCellType>,
491 expected_volume: f64,
492 ) {
493 let mut volume = 0.0;
494 let gdim = grid.geometry_dim();
495 for t in grid.cell_types() {
496 for cell in grid.entity_iter(*t) {
497 let g = cell.geometry();
498 let mut point = vec![0.0; gdim];
499 let mut min_p = vec![10.0; gdim];
500 let mut max_p = vec![-10.0; gdim];
501 for p in g.points() {
502 p.coords(&mut point);
503 for (j, v) in point.iter().enumerate() {
504 if *v < min_p[j] {
505 min_p[j] = *v;
506 }
507 if *v > max_p[j] {
508 max_p[j] = *v;
509 }
510 }
511 }
512 volume += match cell.entity_type() {
513 ReferenceCellType::Interval => {
514 max(&izip!(min_p, max_p).map(|(i, j)| j - i).collect::<Vec<_>>())
515 }
516 ReferenceCellType::Triangle => match gdim {
517 2 => (max_p[0] - min_p[0]) * (max_p[1] - min_p[1]) / 2.0,
518 3 => max(&[
519 (max_p[0] - min_p[0]) * (max_p[1] - min_p[1]) / 2.0,
520 (max_p[0] - min_p[0]) * (max_p[2] - min_p[2]) / 2.0,
521 (max_p[1] - min_p[1]) * (max_p[2] - min_p[2]) / 2.0,
522 ]),
523 _ => {
524 panic!("Unsupported dimension");
525 }
526 },
527 ReferenceCellType::Quadrilateral => match gdim {
528 2 => (max_p[0] - min_p[0]) * (max_p[1] - min_p[1]),
529 3 => max(&[
530 (max_p[0] - min_p[0]) * (max_p[1] - min_p[1]),
531 (max_p[0] - min_p[0]) * (max_p[2] - min_p[2]),
532 (max_p[1] - min_p[1]) * (max_p[2] - min_p[2]),
533 ]),
534 _ => {
535 panic!("Unsupported dimension");
536 }
537 },
538 ReferenceCellType::Tetrahedron => {
539 (max_p[0] - min_p[0]) * (max_p[1] - min_p[1]) * (max_p[2] - min_p[2]) / 6.0
540 }
541 ReferenceCellType::Hexahedron => {
542 (max_p[0] - min_p[0]) * (max_p[1] - min_p[1]) * (max_p[2] - min_p[2])
543 }
544 _ => {
545 panic!("Unsupported cell");
546 }
547 };
548 }
549 }
550 assert_relative_eq!(volume, expected_volume, epsilon = 1e-10);
551 }
552
553 #[test]
554 fn test_unit_interval() {
555 check_volume(&unit_interval::<f64>(1), 1.0);
556 check_volume(&unit_interval::<f64>(2), 1.0);
557 check_volume(&unit_interval::<f64>(4), 1.0);
558 check_volume(&unit_interval::<f64>(7), 1.0);
559 }
560
561 #[test]
562 fn test_unit_square_triangle() {
563 check_volume(&unit_square::<f64>(1, 1, ReferenceCellType::Triangle), 1.0);
564 check_volume(&unit_square::<f64>(2, 2, ReferenceCellType::Triangle), 1.0);
565 check_volume(&unit_square::<f64>(4, 5, ReferenceCellType::Triangle), 1.0);
566 check_volume(&unit_square::<f64>(7, 6, ReferenceCellType::Triangle), 1.0);
567 }
568
569 #[test]
570 fn test_unit_square_quadrilateral() {
571 check_volume(
572 &unit_square::<f64>(1, 1, ReferenceCellType::Quadrilateral),
573 1.0,
574 );
575 check_volume(
576 &unit_square::<f64>(2, 2, ReferenceCellType::Quadrilateral),
577 1.0,
578 );
579 check_volume(
580 &unit_square::<f64>(4, 5, ReferenceCellType::Quadrilateral),
581 1.0,
582 );
583 check_volume(
584 &unit_square::<f64>(7, 6, ReferenceCellType::Quadrilateral),
585 1.0,
586 );
587 }
588
589 #[test]
590 fn test_unit_square_boundary() {
591 check_volume(&unit_square_boundary::<f64>(1, 1), 4.0);
592 check_volume(&unit_square_boundary::<f64>(2, 2), 4.0);
593 check_volume(&unit_square_boundary::<f64>(4, 5), 4.0);
594 check_volume(&unit_square_boundary::<f64>(7, 6), 4.0);
595 }
596
597 #[test]
598 fn test_unit_cube_boundary_triangle() {
599 check_volume(
600 &unit_cube_boundary::<f64>(1, 1, 1, ReferenceCellType::Triangle),
601 6.0,
602 );
603 check_volume(
604 &unit_cube_boundary::<f64>(2, 2, 2, ReferenceCellType::Triangle),
605 6.0,
606 );
607 check_volume(
608 &unit_cube_boundary::<f64>(4, 5, 5, ReferenceCellType::Triangle),
609 6.0,
610 );
611 check_volume(
612 &unit_cube_boundary::<f64>(7, 6, 4, ReferenceCellType::Triangle),
613 6.0,
614 );
615 }
616
617 #[test]
618 fn test_unit_cube_boundary_quadrilateral() {
619 check_volume(
620 &unit_cube_boundary::<f64>(1, 1, 1, ReferenceCellType::Quadrilateral),
621 6.0,
622 );
623 check_volume(
624 &unit_cube_boundary::<f64>(2, 2, 2, ReferenceCellType::Quadrilateral),
625 6.0,
626 );
627 check_volume(
628 &unit_cube_boundary::<f64>(4, 5, 5, ReferenceCellType::Quadrilateral),
629 6.0,
630 );
631 check_volume(
632 &unit_cube_boundary::<f64>(7, 6, 4, ReferenceCellType::Quadrilateral),
633 6.0,
634 );
635 }
636
637 #[test]
638 fn test_unit_cube_tetrahedron() {
639 check_volume(
640 &unit_cube::<f64>(1, 1, 1, ReferenceCellType::Tetrahedron),
641 1.0,
642 );
643 check_volume(
644 &unit_cube::<f64>(2, 2, 2, ReferenceCellType::Tetrahedron),
645 1.0,
646 );
647 check_volume(
648 &unit_cube::<f64>(4, 5, 5, ReferenceCellType::Tetrahedron),
649 1.0,
650 );
651 check_volume(
652 &unit_cube::<f64>(7, 6, 4, ReferenceCellType::Tetrahedron),
653 1.0,
654 );
655 }
656 #[test]
657 fn test_unit_cube_hexahedron() {
658 check_volume(
659 &unit_cube::<f64>(1, 1, 1, ReferenceCellType::Hexahedron),
660 1.0,
661 );
662 check_volume(
663 &unit_cube::<f64>(2, 2, 2, ReferenceCellType::Hexahedron),
664 1.0,
665 );
666 check_volume(
667 &unit_cube::<f64>(4, 5, 5, ReferenceCellType::Hexahedron),
668 1.0,
669 );
670 check_volume(
671 &unit_cube::<f64>(7, 6, 4, ReferenceCellType::Hexahedron),
672 1.0,
673 );
674 }
675
676 #[test]
677 fn test_unit_cube_edges() {
678 check_volume(&unit_cube_edges::<f64>(1, 1, 1), 12.0);
679 check_volume(&unit_cube_edges::<f64>(2, 2, 2), 12.0);
680 check_volume(&unit_cube_edges::<f64>(4, 5, 5), 12.0);
681 check_volume(&unit_cube_edges::<f64>(7, 6, 4), 12.0);
682 }
683}