1use crate::{
2 HasTopology,
3 element::{Handle, VH},
4 error::Error,
5 iterator::HasIterators,
6 mesh::{Adaptor, FloatScalarAdaptor, PolyMeshT},
7};
8use std::ops::{Add, Div, Mul, Neg};
9
10impl<A> PolyMeshT<3, A>
11where
12 A: Adaptor<3>,
13{
14 pub fn quad_box(min: A::Vector, max: A::Vector) -> Result<Self, Error> {
29 const BOX_POS: [(bool, bool, bool); 8] = [
30 (false, false, false),
31 (true, false, false),
32 (true, true, false),
33 (false, true, false),
34 (false, false, true),
35 (true, false, true),
36 (true, true, true),
37 (false, true, true),
38 ];
39 const BOX_IDX: [(u32, u32, u32, u32); 6] = [
40 (0, 3, 2, 1),
41 (0, 1, 5, 4),
42 (1, 2, 6, 5),
43 (2, 3, 7, 6),
44 (3, 0, 4, 7),
45 (4, 5, 6, 7),
46 ];
47 let mut qbox = Self::with_capacity(8, 12, 6);
48 let mut pos: [A::Vector; 8] = [A::zero_vector(); 8];
49 for (i, &(xf, yf, zf)) in BOX_POS.iter().enumerate() {
50 pos[i] = A::vector([
51 A::vector_coord(if xf { &max } else { &min }, 0),
52 A::vector_coord(if yf { &max } else { &min }, 1),
53 A::vector_coord(if zf { &max } else { &min }, 2),
54 ]);
55 }
56 let verts = qbox.add_vertices(&pos)?;
57 assert_eq!(
58 verts,
59 (0..8).into(),
60 "Vertices are expected to be in one contiguous range"
61 );
62 for (a, b, c, d) in BOX_IDX {
63 qbox.add_quad_face(a.into(), b.into(), c.into(), d.into())?;
64 }
65 Ok(qbox)
66 }
67
68 pub fn unit_box() -> Result<Self, Error>
71 where
72 A: FloatScalarAdaptor<3>,
73 {
74 Self::quad_box(
75 A::vector([A::scalarf64(0.0); 3]),
76 A::vector([A::scalarf64(1.0); 3]),
77 )
78 }
79}
80
81impl<const DIM: usize, A> PolyMeshT<DIM, A>
82where
83 A: Adaptor<DIM> + FloatScalarAdaptor<DIM>,
84 A::Scalar: Add<Output = A::Scalar>,
85 A::Vector: Add<Output = A::Vector> + Div<A::Scalar, Output = A::Vector>,
86{
87 pub fn dual_mesh(&self) -> Result<Self, Error> {
98 let mut outmesh =
99 Self::with_capacity(self.num_faces(), self.num_edges(), self.num_vertices());
100 let points = self.points();
101 let points = points.try_borrow()?;
102 let fstatus = self.topol.fstatus.try_borrow()?;
103 for f in self.faces() {
104 if fstatus[f].deleted() {
105 return Err(Error::DeletedFace(f));
106 }
107 outmesh.add_vertex(self.calc_face_centroid(f, &points))?;
108 }
109 let mut fverts: Vec<VH> = Vec::new();
110 for v in self.vertices().filter(|v| !v.is_boundary(self)) {
111 fverts.clear();
112 fverts.extend(self.vf_ccw_iter(v).map(|f| -> VH { f.index().into() }));
113 outmesh.add_face(&fverts)?;
114 }
115 Ok(outmesh)
116 }
117}
118
119impl<A> PolyMeshT<3, A>
121where
122 A: Adaptor<3> + FloatScalarAdaptor<3>,
123 A::Scalar: Mul<Output = A::Scalar> + Neg<Output = A::Scalar>,
124{
125 pub fn tetrahedron(radius: A::Scalar) -> Result<Self, Error> {
128 let mut mesh = Self::with_capacity(4, 6, 4);
129 let a = radius * A::scalarf64(1.0f64 / 3.0);
130 let b = radius * A::scalarf64((8.0 / 9.0f64).sqrt());
131 let c = radius * A::scalarf64((2.0 / 9.0f64).sqrt());
132 let d = radius * A::scalarf64((2.0 / 3.0f64).sqrt());
133 let verts = mesh.add_vertices(&[
134 A::vector([A::scalarf64(0.0), A::scalarf64(0.0), A::scalarf64(1.0)]),
135 A::vector([-c, d, -a]),
136 A::vector([-c, -d, -a]),
137 A::vector([b, A::scalarf64(0.0), -a]),
138 ])?;
139 assert_eq!(
140 verts,
141 (0..4).into(),
142 "Vertices are expected to be in one contiguous range"
143 );
144 mesh.add_tri_face(0.into(), 1.into(), 2.into())?;
145 mesh.add_tri_face(0.into(), 2.into(), 3.into())?;
146 mesh.add_tri_face(0.into(), 3.into(), 1.into())?;
147 mesh.add_tri_face(3.into(), 2.into(), 1.into())?;
148 Ok(mesh)
149 }
150
151 pub fn hexahedron(radius: A::Scalar) -> Result<Self, Error> {
154 let a = radius * A::scalarf64(1.0f64 / 3.0f64.sqrt());
155 let mut mesh = Self::with_capacity(8, 12, 6);
156 let verts = mesh.add_vertices(&[
157 A::vector([-a, -a, -a]),
158 A::vector([a, -a, -a]),
159 A::vector([a, a, -a]),
160 A::vector([-a, a, -a]),
161 A::vector([-a, -a, a]),
162 A::vector([a, -a, a]),
163 A::vector([a, a, a]),
164 A::vector([-a, a, a]),
165 ])?;
166 assert_eq!(
167 verts,
168 (0..8).into(),
169 "Vertices are expected to be in one contiguous range"
170 );
171 mesh.add_quad_face(3.into(), 2.into(), 1.into(), 0.into())?;
172 mesh.add_quad_face(2.into(), 6.into(), 5.into(), 1.into())?;
173 mesh.add_quad_face(5.into(), 6.into(), 7.into(), 4.into())?;
174 mesh.add_quad_face(0.into(), 4.into(), 7.into(), 3.into())?;
175 mesh.add_quad_face(3.into(), 7.into(), 6.into(), 2.into())?;
176 mesh.add_quad_face(1.into(), 5.into(), 4.into(), 0.into())?;
177 Ok(mesh)
178 }
179
180 pub fn octahedron(radius: A::Scalar) -> Result<Self, Error> {
183 let mut mesh = Self::with_capacity(6, 12, 8);
184 let zero = A::scalarf64(0.0);
185 let verts = mesh.add_vertices(&[
186 A::vector([radius, zero, zero]),
187 A::vector([zero, radius, zero]),
188 A::vector([-radius, zero, zero]),
189 A::vector([zero, -radius, zero]),
190 A::vector([zero, zero, radius]),
191 A::vector([zero, zero, -radius]),
192 ])?;
193 assert_eq!(
194 verts,
195 (0..6).into(),
196 "Vertices are expected to be in one contiguous range"
197 );
198 mesh.add_tri_face(0.into(), 4.into(), 3.into())?;
199 mesh.add_tri_face(1.into(), 4.into(), 0.into())?;
200 mesh.add_tri_face(2.into(), 4.into(), 1.into())?;
201 mesh.add_tri_face(3.into(), 4.into(), 2.into())?;
202 mesh.add_tri_face(3.into(), 5.into(), 0.into())?;
203 mesh.add_tri_face(0.into(), 5.into(), 1.into())?;
204 mesh.add_tri_face(1.into(), 5.into(), 2.into())?;
205 mesh.add_tri_face(2.into(), 5.into(), 3.into())?;
206 Ok(mesh)
207 }
208
209 pub fn icosahedron(radius: A::Scalar) -> Result<Self, Error> {
212 let mut mesh = Self::with_capacity(12, 30, 20);
213 let verts = mesh.add_vertices(&[
214 A::vector([
215 radius * A::scalarf64(0.0),
216 radius * A::scalarf64(0.5257311121191336),
217 radius * A::scalarf64(-0.8506508083520399),
218 ]),
219 A::vector([
220 radius * A::scalarf64(0.5257311121191336),
221 radius * A::scalarf64(0.8506508083520399),
222 radius * A::scalarf64(0.0),
223 ]),
224 A::vector([
225 radius * A::scalarf64(-0.5257311121191336),
226 radius * A::scalarf64(0.8506508083520399),
227 radius * A::scalarf64(0.0),
228 ]),
229 A::vector([
230 radius * A::scalarf64(0.0),
231 radius * A::scalarf64(0.5257311121191336),
232 radius * A::scalarf64(0.8506508083520399),
233 ]),
234 A::vector([
235 radius * A::scalarf64(0.0),
236 radius * A::scalarf64(-0.5257311121191336),
237 radius * A::scalarf64(0.8506508083520399),
238 ]),
239 A::vector([
240 radius * A::scalarf64(-0.8506508083520399),
241 radius * A::scalarf64(0.0),
242 radius * A::scalarf64(0.5257311121191336),
243 ]),
244 A::vector([
245 radius * A::scalarf64(0.0),
246 radius * A::scalarf64(-0.5257311121191336),
247 radius * A::scalarf64(-0.8506508083520399),
248 ]),
249 A::vector([
250 radius * A::scalarf64(0.8506508083520399),
251 radius * A::scalarf64(0.0),
252 radius * A::scalarf64(-0.5257311121191336),
253 ]),
254 A::vector([
255 radius * A::scalarf64(0.8506508083520399),
256 radius * A::scalarf64(0.0),
257 radius * A::scalarf64(0.5257311121191336),
258 ]),
259 A::vector([
260 radius * A::scalarf64(-0.8506508083520399),
261 radius * A::scalarf64(0.0),
262 radius * A::scalarf64(-0.5257311121191336),
263 ]),
264 A::vector([
265 radius * A::scalarf64(0.5257311121191336),
266 radius * A::scalarf64(-0.8506508083520399),
267 radius * A::scalarf64(0.0),
268 ]),
269 A::vector([
270 radius * A::scalarf64(-0.5257311121191336),
271 radius * A::scalarf64(-0.8506508083520399),
272 radius * A::scalarf64(0.0),
273 ]),
274 ])?;
275 assert_eq!(
276 verts,
277 (0..12).into(),
278 "Vertices are expected to be in one contiguous range"
279 );
280 mesh.add_tri_face(2.into(), 1.into(), 0.into())?;
281 mesh.add_tri_face(1.into(), 2.into(), 3.into())?;
282 mesh.add_tri_face(5.into(), 4.into(), 3.into())?;
283 mesh.add_tri_face(4.into(), 8.into(), 3.into())?;
284 mesh.add_tri_face(7.into(), 6.into(), 0.into())?;
285 mesh.add_tri_face(6.into(), 9.into(), 0.into())?;
286 mesh.add_tri_face(11.into(), 10.into(), 4.into())?;
287 mesh.add_tri_face(10.into(), 11.into(), 6.into())?;
288 mesh.add_tri_face(9.into(), 5.into(), 2.into())?;
289 mesh.add_tri_face(5.into(), 9.into(), 11.into())?;
290 mesh.add_tri_face(8.into(), 7.into(), 1.into())?;
291 mesh.add_tri_face(7.into(), 8.into(), 10.into())?;
292 mesh.add_tri_face(2.into(), 5.into(), 3.into())?;
293 mesh.add_tri_face(8.into(), 1.into(), 3.into())?;
294 mesh.add_tri_face(9.into(), 2.into(), 0.into())?;
295 mesh.add_tri_face(1.into(), 7.into(), 0.into())?;
296 mesh.add_tri_face(11.into(), 9.into(), 6.into())?;
297 mesh.add_tri_face(7.into(), 10.into(), 6.into())?;
298 mesh.add_tri_face(5.into(), 11.into(), 4.into())?;
299 mesh.add_tri_face(10.into(), 8.into(), 4.into())?;
300 Ok(mesh)
301 }
302
303 pub fn dodecahedron(radius: A::Scalar) -> Result<Self, Error> {
306 let mut mesh = Self::with_capacity(20, 30, 12);
307 let verts = mesh.add_vertices(&[
308 A::vector([
309 radius * A::scalarf64(0.0),
310 radius * A::scalarf64(0.9341723589627157),
311 radius * A::scalarf64(-0.35682208977308993),
312 ]),
313 A::vector([
314 radius * A::scalarf64(0.0),
315 radius * A::scalarf64(0.9341723589627157),
316 radius * A::scalarf64(0.35682208977308993),
317 ]),
318 A::vector([
319 radius * A::scalarf64(-0.35682208977308993),
320 radius * A::scalarf64(0.0),
321 radius * A::scalarf64(0.9341723589627157),
322 ]),
323 A::vector([
324 radius * A::scalarf64(0.35682208977308993),
325 radius * A::scalarf64(0.0),
326 radius * A::scalarf64(0.9341723589627157),
327 ]),
328 A::vector([
329 radius * A::scalarf64(0.35682208977308993),
330 radius * A::scalarf64(0.0),
331 radius * A::scalarf64(-0.9341723589627157),
332 ]),
333 A::vector([
334 radius * A::scalarf64(-0.35682208977308993),
335 radius * A::scalarf64(0.0),
336 radius * A::scalarf64(-0.9341723589627157),
337 ]),
338 A::vector([
339 radius * A::scalarf64(0.0),
340 radius * A::scalarf64(-0.9341723589627157),
341 radius * A::scalarf64(0.35682208977308993),
342 ]),
343 A::vector([
344 radius * A::scalarf64(0.0),
345 radius * A::scalarf64(-0.9341723589627157),
346 radius * A::scalarf64(-0.35682208977308993),
347 ]),
348 A::vector([
349 radius * A::scalarf64(-0.9341723589627157),
350 radius * A::scalarf64(0.35682208977308993),
351 radius * A::scalarf64(0.0),
352 ]),
353 A::vector([
354 radius * A::scalarf64(-0.9341723589627157),
355 radius * A::scalarf64(-0.35682208977308993),
356 radius * A::scalarf64(0.0),
357 ]),
358 A::vector([
359 radius * A::scalarf64(0.9341723589627157),
360 radius * A::scalarf64(0.35682208977308993),
361 radius * A::scalarf64(0.0),
362 ]),
363 A::vector([
364 radius * A::scalarf64(0.9341723589627157),
365 radius * A::scalarf64(-0.35682208977308993),
366 radius * A::scalarf64(0.0),
367 ]),
368 A::vector([
369 radius * A::scalarf64(-0.5773502691896257),
370 radius * A::scalarf64(0.5773502691896257),
371 radius * A::scalarf64(0.5773502691896257),
372 ]),
373 A::vector([
374 radius * A::scalarf64(0.5773502691896257),
375 radius * A::scalarf64(0.5773502691896257),
376 radius * A::scalarf64(0.5773502691896257),
377 ]),
378 A::vector([
379 radius * A::scalarf64(-0.5773502691896257),
380 radius * A::scalarf64(0.5773502691896257),
381 radius * A::scalarf64(-0.5773502691896257),
382 ]),
383 A::vector([
384 radius * A::scalarf64(0.5773502691896257),
385 radius * A::scalarf64(0.5773502691896257),
386 radius * A::scalarf64(-0.5773502691896257),
387 ]),
388 A::vector([
389 radius * A::scalarf64(-0.5773502691896257),
390 radius * A::scalarf64(-0.5773502691896257),
391 radius * A::scalarf64(-0.5773502691896257),
392 ]),
393 A::vector([
394 radius * A::scalarf64(0.5773502691896257),
395 radius * A::scalarf64(-0.5773502691896257),
396 radius * A::scalarf64(-0.5773502691896257),
397 ]),
398 A::vector([
399 radius * A::scalarf64(-0.5773502691896257),
400 radius * A::scalarf64(-0.5773502691896257),
401 radius * A::scalarf64(0.5773502691896257),
402 ]),
403 A::vector([
404 radius * A::scalarf64(0.5773502691896257),
405 radius * A::scalarf64(-0.5773502691896257),
406 radius * A::scalarf64(0.5773502691896257),
407 ]),
408 ])?;
409 assert_eq!(
410 verts,
411 (0..20).into(),
412 "Vertices are expected to be in one contiguous range"
413 );
414 mesh.add_face(&[15.into(), 4.into(), 5.into(), 14.into(), 0.into()])?;
415 mesh.add_face(&[15.into(), 0.into(), 1.into(), 13.into(), 10.into()])?;
416 mesh.add_face(&[14.into(), 8.into(), 12.into(), 1.into(), 0.into()])?;
417 mesh.add_face(&[13.into(), 1.into(), 12.into(), 2.into(), 3.into()])?;
418 mesh.add_face(&[19.into(), 3.into(), 2.into(), 18.into(), 6.into()])?;
419 mesh.add_face(&[18.into(), 2.into(), 12.into(), 8.into(), 9.into()])?;
420 mesh.add_face(&[17.into(), 7.into(), 16.into(), 5.into(), 4.into()])?;
421 mesh.add_face(&[17.into(), 4.into(), 15.into(), 10.into(), 11.into()])?;
422 mesh.add_face(&[19.into(), 11.into(), 10.into(), 13.into(), 3.into()])?;
423 mesh.add_face(&[16.into(), 9.into(), 8.into(), 14.into(), 5.into()])?;
424 mesh.add_face(&[19.into(), 6.into(), 7.into(), 17.into(), 11.into()])?;
425 mesh.add_face(&[18.into(), 9.into(), 16.into(), 7.into(), 6.into()])?;
426 Ok(mesh)
427 }
428}
429
430#[cfg(test)]
431mod test {
432 use crate::{
433 HasTopology,
434 iterator::HasIterators,
435 macros::assert_f32_eq,
436 mesh::{PolyMeshF32, Vec3},
437 };
438 use core::f32;
439
440 #[test]
441 fn t_quad_box() {
442 let qbox = PolyMeshF32::quad_box(Vec3(0., 0., 0.), Vec3(1., 1., 1.))
443 .expect("Cannot create a quad box mesh");
444 assert_eq!(qbox.num_vertices(), 8);
445 assert_eq!(qbox.num_halfedges(), 24);
446 assert_eq!(qbox.num_edges(), 12);
447 assert_eq!(qbox.num_faces(), 6);
448 for v in qbox.vertices() {
449 assert_eq!(qbox.vf_ccw_iter(v).count(), 3);
450 }
451 assert_eq!(1., qbox.try_calc_volume().expect("Cannot compute volume"));
452 assert_eq!(6., qbox.try_calc_area().expect("Cannot compute area"));
453 }
454
455 #[test]
456 fn t_tetrahedron() {
457 let tet = PolyMeshF32::tetrahedron(1.0).expect("Cannot create a tetrahedron");
458 assert_eq!(4, tet.num_vertices());
459 assert_eq!(12, tet.num_halfedges());
460 assert_eq!(6, tet.num_edges());
461 assert_eq!(4, tet.num_faces());
462 assert_eq!(
463 8.0 / 3.0f32.sqrt(),
464 tet.try_calc_area().expect("Cannot compute area")
465 );
466 assert_f32_eq!(
467 8.0 / (9.0 * 3.0f32.sqrt()),
468 tet.try_calc_volume().expect("Cannot compute volume")
469 );
470 }
471
472 #[test]
473 fn t_hexahedron() {
474 let hex = PolyMeshF32::hexahedron(1.0).expect("Cannot create hexahedron");
475 assert_eq!(hex.num_vertices(), 8);
476 assert_eq!(hex.num_halfedges(), 24);
477 assert_eq!(hex.num_edges(), 12);
478 assert_eq!(hex.num_faces(), 6);
479 assert_f32_eq!(8.0, hex.try_calc_area().expect("Cannot compute area"), 1e-6);
480 assert_f32_eq!(
481 8.0 / (3.0 * 3.0f32.sqrt()),
482 hex.try_calc_volume().expect("Cannot compute volume")
483 );
484 }
485
486 #[test]
487 fn t_octahedron() {
488 let oct = PolyMeshF32::octahedron(1.0).expect("Cannot create octahedron");
489 assert_eq!(oct.num_vertices(), 6);
490 assert_eq!(oct.num_halfedges(), 24);
491 assert_eq!(oct.num_edges(), 12);
492 assert_eq!(oct.num_faces(), 8);
493 assert_eq!(
494 4.0 * 3.0f32.sqrt(),
495 oct.try_calc_area().expect("Cannot compute area")
496 );
497 assert_f32_eq!(
498 4.0 / 3.0,
499 oct.try_calc_volume().expect("Cannot compute volume")
500 );
501 }
502
503 #[test]
504 fn t_icosahedron() {
505 let ico = PolyMeshF32::icosahedron(1.0).expect("Cannot create icosahedron");
506 assert_eq!(12, ico.num_vertices());
507 assert_eq!(60, ico.num_halfedges());
508 assert_eq!(30, ico.num_edges());
509 assert_eq!(20, ico.num_faces());
510 assert_f32_eq!(
511 {
512 let phi = (1.0 + 5.0f32.sqrt()) / 2.0;
513 20.0 * 3.0f32.sqrt() / (phi * phi + 1.0)
514 },
515 ico.try_calc_area().expect("Cannot compute area"),
516 1e-6
517 );
518 assert_f32_eq!(
519 {
520 let phi = (1.0 + 5.0f32.sqrt()) / 2.0;
521 20.0 * phi * phi / (3.0 * (phi * phi + 1.0) * (phi * phi + 1.0).sqrt())
522 },
523 ico.try_calc_volume().expect("Cannot compute volume"),
524 1e-6
525 );
526 }
527
528 #[test]
529 fn t_dodecahedron() {
530 let dod = PolyMeshF32::dodecahedron(1.0).expect("Cannot create dodecahedron");
531 assert_eq!(20, dod.num_vertices());
532 assert_eq!(60, dod.num_halfedges());
533 assert_eq!(30, dod.num_edges());
534 assert_eq!(12, dod.num_faces());
535 assert_f32_eq!(
536 {
537 let phi = (1.0 + 5.0f32.sqrt()) / 2.0;
538 20.0f32 / (phi * (3.0f32 - phi).sqrt())
539 },
540 dod.try_calc_area().expect("Cannot compute area"),
541 1e-6
542 );
543 assert_f32_eq!(
544 {
545 let phi = (1.0 + 5.0f32.sqrt()) / 2.0;
546 40.0 / (3.0 * 3.0f32.sqrt() * (6.0 - 2.0 * phi))
547 },
548 dod.try_calc_volume().expect("Cannot compute volume")
549 );
550 }
551}