alum/
create.rs

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    /// Makes a box with the following topology, spanning from the min point to
15    /// the max point.
16    ///
17    ///  ```text
18    ///       7-----------6
19    ///      /|          /|
20    ///     / |         / |
21    ///    4-----------5  |
22    ///    |  |        |  |
23    ///    |  3--------|--2
24    ///    | /         | /
25    ///    |/          |/
26    ///    0-----------1
27    ///  ```
28    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    /// Create a mesh representing a box with quadrilateral faces, of size 1,
69    /// spanning from the origin to (1, 1, 1).
70    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    /// Create a mesh that is the dual of this mesh.
88    ///
89    /// <https://en.wikipedia.org/wiki/Dual_polyhedron>.
90    ///
91    /// The dual mesh will contain vertices at the centroids of the input mesh,
92    /// and an edge connecting the vertices for every pair of adjacent faces in
93    /// the input mesh. The output mesh will contain as many vertices as the
94    /// faces of the input mesh, as many edges as the interior edges of the
95    /// input mesh, and as many faces as the interior vertices of the input
96    /// mesh.
97    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
119/// Platonic solids.
120impl<A> PolyMeshT<3, A>
121where
122    A: Adaptor<3> + FloatScalarAdaptor<3>,
123    A::Scalar: Mul<Output = A::Scalar> + Neg<Output = A::Scalar>,
124{
125    /// Create a tetrahedron centered at the origin, with a unit
126    /// circumradius. The vertices of the mesh will lie on the unit sphere.
127    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    /// Create a hexehedron centered at the origin, with a unit
152    /// circumradius. The vertices of the mesh will lie on the unit sphere.
153    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    /// Create an octahedron centered at the origin, with a unit
181    /// circumradius. The vertices of the mesh will lie on the unit sphere.
182    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    /// Create an icosahedron centered at the origin, with a unit
210    /// circumradius. The vertices of the mesh will lie on the unit sphere.
211    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    /// Create a dodecahedron centered at the origin, with a unit
304    /// circumradius. The vertices of the mesh will lie on the unit sphere.
305    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}