gemlab/mesh/
convert_2d.rs

1use super::{Cell, Features, Mesh, Point};
2use crate::prelude::PointId;
3use crate::shapes::{GeoClass, GeoKind, Scratchpad};
4use crate::StrError;
5use russell_lab::Vector;
6use std::collections::HashMap;
7
8impl Mesh {
9    /// Upgrades or downgrades a mesh with triangles or quadrilaterals
10    ///
11    /// # Notes
12    ///
13    /// 1. All cells must have the same GeoKind
14    /// 2. Only [GeoClass::Tri] and [GeoClass::Qua] are allowed
15    /// 3. The points will be completely renumbered
16    /// 4. The corner tags will be replicated into the new mesh
17    /// 5. The points at the middle of edges will inherit the tag of ONE corresponding middle point
18    pub fn convert_2d(&self, target: GeoKind) -> Result<Mesh, StrError> {
19        //        2,
20        //  s     | ',
21        //        |   ',
22        //  i    10     9,
23        //        |       ',
24        //  d     |         ',   side = 1
25        //        5    14     4,
26        //  e     |             ',
27        //        |               ',         NOTE that two neighboring triangles
28        //  =    11    12    13     8,       will have reversed pairs of local
29        //        |                   ',     points, e.g. (6,7) <-> (7,6)
30        //  2     |      side = 0       ',   Thus, the first point along the
31        //        0-----6-----3-----7-----1  edge of one triangle will be the
32        //        1-----7-----3-----6-----0  second point along the edge of the
33        //  2     |      side = 0       ,'   neighbor triangle, and vice-versa
34        //        |                   ,'     Solution: Flip Normals
35        //  =     8    13    12    11'
36        //        |               ,'
37        //  e     |             ,'
38        //        4    14     5'
39        //  d     |         ,'   side = 1
40        //        |       ,'
41        //  i     9    10'
42        //        |   ,'
43        //  s     | ,'
44        //        2'
45        if self.ndim != 2 {
46            return Err("ndim must be equal to 2");
47        }
48
49        // current number of cells
50        let ncell = self.cells.len();
51        if ncell < 1 {
52            return Err("the conversion requires at least one cell");
53        }
54
55        // get info about first cell in the mesh
56        let source = self.cells[0].kind;
57        let class = source.class();
58        let source_nnode = source.nnode();
59        let target_nnode = target.nnode();
60        if target.class() != class {
61            return Err("target class must equal the GeoClass of current cells");
62        }
63        if target.class() != GeoClass::Tri && target.class() != GeoClass::Qua {
64            return Err("target GeoClass must be Tri or Qua");
65        }
66
67        // info about the edges and corners
68        let source_edge_nnode = source.edge_nnode();
69        let target_edge_nnode = target.edge_nnode();
70        let nedge = target.nedge();
71        let ncorner = if class == GeoClass::Tri { 3 } else { 4 };
72        let mut edge_point_markers = vec![0; nedge];
73
74        // maps target local node id to edge number (to set markers of points on edges)
75        const INTERIOR_OR_CORNER: usize = usize::MAX;
76        let mut target_node_to_edge = vec![INTERIOR_OR_CORNER; target_nnode];
77        for e in 0..nedge {
78            for i in 2..target_edge_nnode {
79                let m = target.edge_node_id(e, i);
80                target_node_to_edge[m] = e;
81            }
82        }
83
84        // to search neighbors in the original mesh
85        let source_features = Features::new(self, true);
86
87        // flag indicating that a connectivity point has not been set
88        const UNSET: usize = usize::MAX;
89        assert!(ncell * target_nnode < UNSET);
90
91        // zeroed new cell
92        let zero_cell = Cell {
93            id: 0,
94            marker: 0,
95            kind: target,
96            points: vec![UNSET; target_nnode],
97        };
98
99        // allocate destination mesh
100        let mut dest = Mesh {
101            ndim: self.ndim,
102            points: Vec::new(),
103            cells: vec![zero_cell; ncell],
104            marked_edges: Vec::new(),
105            marked_faces: Vec::new(),
106        };
107
108        // scratchpad for point interpolation (based on the original mesh)
109        let mut pad = Scratchpad::new(self.ndim, source).unwrap();
110
111        // coordinates of new points
112        let mut x = Vector::new(self.ndim);
113
114        // maps old point id to new point id
115        let mut corners: HashMap<PointId, PointId> = HashMap::new();
116
117        // upgrade or downgrade mesh
118        for cell_id in 0..ncell {
119            // check GeoKind
120            if self.cells[cell_id].kind != source {
121                return Err("all cells must have the same GeoKind");
122            }
123
124            // set the coordinates matrix of the original mesh
125            for m in 0..source_nnode {
126                let p = self.cells[cell_id].points[m];
127                for j in 0..self.ndim {
128                    pad.set_xx(m, j, self.points[p].coords[j]);
129                }
130            }
131
132            // find markers of points on edges
133            if source_edge_nnode > 2 {
134                for e in 0..nedge {
135                    for i in 2..source_edge_nnode {
136                        let m = source.edge_node_id(e, i);
137                        let p = self.cells[cell_id].points[m];
138                        edge_point_markers[e] = self.points[p].marker;
139                    }
140                }
141            }
142
143            // set the new cell data
144            dest.cells[cell_id].id = cell_id;
145            dest.cells[cell_id].marker = self.cells[cell_id].marker;
146
147            // handle corner nodes
148            for m in 0..ncorner {
149                let old_point_id = self.cells[cell_id].points[m];
150                if let Some(new_point_id) = corners.get(&old_point_id) {
151                    dest.cells[cell_id].points[m] = *new_point_id;
152                } else {
153                    let new_point_id = dest.points.len();
154                    dest.points.push(Point {
155                        id: new_point_id,
156                        marker: self.points[old_point_id].marker,
157                        coords: self.points[old_point_id].coords.clone(),
158                    });
159                    dest.cells[cell_id].points[m] = new_point_id;
160                    corners.insert(old_point_id, new_point_id);
161                }
162            }
163
164            // consult neighbors to see if there are points (in the middle of the edge) set already
165            if target_edge_nnode > 2 {
166                let neighbors = source_features.get_neighbors_2d(cell_id);
167                for (e, neigh_cell_id, neigh_e) in neighbors {
168                    // only deal with the centre edge points, not the corner ones (start at 2)
169                    for i in 2..target_edge_nnode {
170                        let n = target.edge_node_id(neigh_e, i); // outward normal
171                        let p = dest.cells[neigh_cell_id].points[n];
172                        if p != UNSET {
173                            let m = target.edge_node_id_inward(e, i); // inward normal
174                            dest.cells[cell_id].points[m] = p;
175                        }
176                    }
177                }
178            }
179
180            // add new points (skip corner nodes)
181            for m in ncorner..target_nnode {
182                if dest.cells[cell_id].points[m] == UNSET {
183                    let e = target_node_to_edge[m];
184                    let marker = if e == INTERIOR_OR_CORNER {
185                        0
186                    } else {
187                        edge_point_markers[e]
188                    };
189                    let new_point_id = dest.points.len();
190                    pad.calc_coords(&mut x, target.reference_coords(m)).unwrap();
191                    dest.points.push(Point {
192                        id: new_point_id,
193                        marker,
194                        coords: x.as_data().clone(),
195                    });
196                    dest.cells[cell_id].points[m] = new_point_id;
197                }
198            }
199        }
200        Ok(dest)
201    }
202}
203
204////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
205
206#[cfg(test)]
207mod tests {
208    use crate::mesh::{Cell, Draw, Mesh, Point, Samples};
209    use crate::shapes::GeoKind;
210    use russell_lab::array_approx_eq;
211
212    const SAVE_FIGURE: bool = false;
213
214    fn draw(mesh: &Mesh, larger: bool, filename: &str) {
215        let mut draw = Draw::new();
216        draw.show_cell_ids(true).show_point_ids(true);
217        if larger {
218            draw.set_size(600.0, 600.0);
219        }
220        draw.all(&mesh, filename).unwrap();
221    }
222
223    #[test]
224    fn convert_mesh_2d_captures_errors() {
225        #[rustfmt::skip]
226        let mesh = Mesh {
227            ndim: 3,
228            points: vec![
229                Point { id: 0, marker: 0, coords: vec![0.0, 0.0, 0.0] },
230                Point { id: 1, marker: 0, coords: vec![1.0, 0.0, 0.0] },
231                Point { id: 2, marker: 0, coords: vec![1.0, 1.0, 0.0] },
232                Point { id: 3, marker: 0, coords: vec![0.0, 1.0, 0.0] },
233                Point { id: 4, marker: 0, coords: vec![0.0, 0.0, 1.0] },
234                Point { id: 5, marker: 0, coords: vec![1.0, 0.0, 1.0] },
235                Point { id: 6, marker: 0, coords: vec![1.0, 1.0, 1.0] },
236                Point { id: 7, marker: 0, coords: vec![0.0, 1.0, 1.0] },
237            ],
238            cells: vec![
239                Cell { id: 0, marker: 1, kind: GeoKind::Hex8, points: vec![0,1,2,3, 4,5,6,7] },
240            ],
241            marked_edges: Vec::new(),
242            marked_faces: Vec::new(),
243        };
244        assert_eq!(mesh.convert_2d(GeoKind::Tri15).err(), Some("ndim must be equal to 2"));
245
246        let mesh = Mesh {
247            ndim: 2,
248            points: vec![],
249            cells: vec![],
250            marked_edges: Vec::new(),
251            marked_faces: Vec::new(),
252        };
253        assert_eq!(
254            mesh.convert_2d(GeoKind::Tri15).err(),
255            Some("the conversion requires at least one cell")
256        );
257
258        #[rustfmt::skip]
259        let mesh = Mesh {
260            ndim: 2,
261            points: vec![
262                Point { id: 0, marker: 0, coords: vec![0.0, 0.0 ] },
263                Point { id: 1, marker: 0, coords: vec![1.0, 0.0 ] },
264                Point { id: 2, marker: 0, coords: vec![0.5, 0.85] },
265            ],
266            cells: vec![
267                Cell { id: 0, marker: 1, kind: GeoKind::Tri3, points: vec![0, 1, 2] },
268            ],
269            marked_edges: Vec::new(),
270            marked_faces: Vec::new(),
271        };
272        assert_eq!(
273            mesh.convert_2d(GeoKind::Qua8).err(),
274            Some("target class must equal the GeoClass of current cells")
275        );
276
277        #[rustfmt::skip]
278        let mesh =Mesh {
279            ndim: 2,
280            points: vec![
281                Point { id: 0, marker: 0, coords: vec![0.0, 0.0] },
282                Point { id: 1, marker: 0, coords: vec![1.0, 1.0] },
283            ],
284            cells: vec![
285                Cell { id: 0, marker: 1, kind: GeoKind::Lin2, points: vec![0, 1] },
286            ],
287            marked_edges: Vec::new(),
288            marked_faces: Vec::new(),
289        };
290        assert_eq!(
291            mesh.convert_2d(GeoKind::Lin3).err(),
292            Some("target GeoClass must be Tri or Qua")
293        );
294
295        #[rustfmt::skip]
296        let mesh = Mesh {
297            ndim: 2,
298            points: vec![
299                Point { id: 0, marker: 0, coords: vec![0.0, 0.0 ] },
300                Point { id: 1, marker: 0, coords: vec![1.0, 0.0 ] },
301                Point { id: 2, marker: 0, coords: vec![0.5, 0.85] },
302                Point { id: 3, marker: 0, coords: vec![1.0, 1.0] },
303            ],
304            cells: vec![
305                Cell { id: 0, marker: 1, kind: GeoKind::Tri3, points: vec![0, 1, 2] },
306                Cell { id: 1, marker: 1, kind: GeoKind::Qua4, points: vec![0, 1, 2, 3] },
307            ],
308            marked_edges: Vec::new(),
309            marked_faces: Vec::new(),
310        };
311        assert_eq!(
312            mesh.convert_2d(GeoKind::Tri6).err(),
313            Some("all cells must have the same GeoKind")
314        );
315    }
316
317    #[test]
318    fn convert_tri6_to_tri15_works() {
319        #[rustfmt::skip]
320        let mesh = Mesh {
321            ndim: 2,
322            points: vec![
323                Point { id: 0, marker: 0, coords:   vec![0.0, 0.0] },
324                Point { id: 1, marker: 0, coords:   vec![4.0, 0.0] },
325                Point { id: 2, marker: 0, coords:   vec![0.0, 4.0] },
326                Point { id: 3, marker: 0, coords:   vec![2.0, 0.0] },
327                Point { id: 4, marker: -20, coords: vec![2.0, 2.0] },
328                Point { id: 5, marker: -10, coords: vec![0.0, 2.0] },
329                Point { id: 6, marker: 0, coords:   vec![0.0,-4.0] },
330                Point { id: 7, marker: -30, coords: vec![2.0,-2.0] },
331                Point { id: 8, marker: -10, coords: vec![0.0,-2.0] },
332            ],
333            cells: vec![
334                Cell { id: 0, marker: 1, kind: GeoKind::Tri6, points: vec![1, 2, 0, 4, 5, 3] },
335                Cell { id: 1, marker: 2, kind: GeoKind::Tri6, points: vec![0, 6, 1, 8, 7, 3] },
336            ],
337            marked_edges: Vec::new(),
338            marked_faces: Vec::new(),
339        };
340        mesh.check_all().unwrap();
341
342        if SAVE_FIGURE {
343            draw(&mesh, false, "/tmp/gemlab/test_convert_2d_tri6_to_tri15_before.svg");
344        }
345
346        let res = mesh.convert_2d(GeoKind::Tri15).unwrap();
347
348        if SAVE_FIGURE {
349            draw(&mesh, false, "/tmp/gemlab/test_convert_2d_tri6_to_tri15_after.svg");
350        }
351
352        res.check_all().unwrap();
353        res.check_overlapping_points(0.02).unwrap();
354
355        assert_eq!(res.points.len(), 2 * 15 - 5);
356        assert_eq!(res.cells[0].points, (0..15).collect::<Vec<_>>());
357        assert_eq!(
358            res.cells[1].points,
359            &[2, 15, 0, 16, 17, 5, 18, 19, 20, 21, 11, 10, 22, 23, 24],
360        );
361
362        array_approx_eq(&res.points[0].coords, &[4.0, 0.0], 1e-15);
363        array_approx_eq(&res.points[1].coords, &[0.0, 4.0], 1e-15);
364        array_approx_eq(&res.points[2].coords, &[0.0, 0.0], 1e-15);
365        array_approx_eq(&res.points[3].coords, &[2.0, 2.0], 1e-15);
366        array_approx_eq(&res.points[4].coords, &[0.0, 2.0], 1e-15);
367        array_approx_eq(&res.points[5].coords, &[2.0, 0.0], 1e-15);
368        array_approx_eq(&res.points[6].coords, &[3.0, 1.0], 1e-15);
369        array_approx_eq(&res.points[7].coords, &[1.0, 3.0], 1e-15);
370        array_approx_eq(&res.points[8].coords, &[0.0, 3.0], 1e-15);
371        array_approx_eq(&res.points[9].coords, &[0.0, 1.0], 1e-15);
372        array_approx_eq(&res.points[10].coords, &[1.0, 0.0], 1e-15);
373        array_approx_eq(&res.points[11].coords, &[3.0, 0.0], 1e-15);
374        array_approx_eq(&res.points[12].coords, &[2.0, 1.0], 1e-15);
375        array_approx_eq(&res.points[13].coords, &[1.0, 2.0], 1e-15);
376        array_approx_eq(&res.points[14].coords, &[1.0, 1.0], 1e-15);
377        array_approx_eq(&res.points[15].coords, &[0.0, -4.0], 1e-15);
378        array_approx_eq(&res.points[16].coords, &[0.0, -2.0], 1e-15);
379        array_approx_eq(&res.points[17].coords, &[2.0, -2.0], 1e-15);
380        array_approx_eq(&res.points[18].coords, &[0.0, -1.0], 1e-15);
381        array_approx_eq(&res.points[19].coords, &[0.0, -3.0], 1e-15);
382        array_approx_eq(&res.points[20].coords, &[1.0, -3.0], 1e-15);
383        array_approx_eq(&res.points[21].coords, &[3.0, -1.0], 1e-15);
384        array_approx_eq(&res.points[22].coords, &[1.0, -1.0], 1e-15);
385        array_approx_eq(&res.points[23].coords, &[1.0, -2.0], 1e-15);
386        array_approx_eq(&res.points[24].coords, &[2.0, -1.0], 1e-15);
387
388        assert_eq!(res.points[0].marker, 0);
389        assert_eq!(res.points[1].marker, 0);
390        assert_eq!(res.points[2].marker, 0);
391        assert_eq!(res.points[3].marker, -20);
392        assert_eq!(res.points[4].marker, -10);
393        assert_eq!(res.points[5].marker, 0);
394        assert_eq!(res.points[6].marker, -20);
395        assert_eq!(res.points[7].marker, -20);
396        assert_eq!(res.points[8].marker, -10);
397        assert_eq!(res.points[9].marker, -10);
398        assert_eq!(res.points[10].marker, 0);
399        assert_eq!(res.points[11].marker, 0);
400        assert_eq!(res.points[12].marker, 0);
401        assert_eq!(res.points[13].marker, 0);
402        assert_eq!(res.points[14].marker, 0);
403        assert_eq!(res.points[15].marker, 0);
404        assert_eq!(res.points[16].marker, -10);
405        assert_eq!(res.points[17].marker, -30);
406        assert_eq!(res.points[18].marker, -10);
407        assert_eq!(res.points[19].marker, -10);
408        assert_eq!(res.points[20].marker, -30);
409        assert_eq!(res.points[21].marker, -30);
410        assert_eq!(res.points[22].marker, 0);
411        assert_eq!(res.points[23].marker, 0);
412        assert_eq!(res.points[24].marker, 0);
413    }
414
415    #[test]
416    fn convert_tri6_to_tri10_works() {
417        #[rustfmt::skip]
418        let mesh = Mesh {
419            ndim: 2,
420            points: vec![
421                Point { id: 0, marker: 0, coords: vec![0.0, 0.0] },
422                Point { id: 1, marker: 0, coords: vec![4.0, 0.0] },
423                Point { id: 2, marker: 0, coords: vec![0.0, 4.0] },
424                Point { id: 3, marker: 0, coords: vec![2.0, 0.0] },
425                Point { id: 4, marker: -20, coords: vec![2.0, 2.0] },
426                Point { id: 5, marker: -10, coords: vec![0.0, 2.0] },
427                Point { id: 6, marker: 0, coords: vec![0.0,-4.0] },
428                Point { id: 7, marker: -30, coords: vec![2.0,-2.0] },
429                Point { id: 8, marker: -10, coords: vec![0.0,-2.0] },
430            ],
431            cells: vec![
432                Cell { id: 0, marker: 1, kind: GeoKind::Tri6, points: vec![1, 2, 0, 4, 5, 3] },
433                Cell { id: 1, marker: 2, kind: GeoKind::Tri6, points: vec![0, 6, 1, 8, 7, 3] },
434            ],
435            marked_edges: Vec::new(),
436            marked_faces: Vec::new(),
437        };
438        mesh.check_all().unwrap();
439
440        if SAVE_FIGURE {
441            draw(&mesh, false, "/tmp/gemlab/test_convert_2d_tri6_to_tri10_before.svg");
442        }
443
444        let res = mesh.convert_2d(GeoKind::Tri10).unwrap();
445
446        if SAVE_FIGURE {
447            draw(&res, false, "/tmp/gemlab/test_convert_2d_tri6_to_tri10_after.svg");
448        }
449
450        res.check_all().unwrap();
451        res.check_overlapping_points(0.02).unwrap();
452
453        assert_eq!(res.points.len(), 16);
454
455        let ft = 4.0 / 3.0;
456        let et = 8.0 / 3.0;
457        array_approx_eq(&res.points[0].coords, &[4.0, 0.0], 1e-15);
458        array_approx_eq(&res.points[1].coords, &[0.0, 4.0], 1e-15);
459        array_approx_eq(&res.points[2].coords, &[0.0, 0.0], 1e-15);
460        array_approx_eq(&res.points[3].coords, &[et, ft], 1e-15);
461        array_approx_eq(&res.points[4].coords, &[0.0, et], 1e-15);
462        array_approx_eq(&res.points[5].coords, &[ft, 0.0], 1e-15);
463        array_approx_eq(&res.points[6].coords, &[ft, et], 1e-15);
464        array_approx_eq(&res.points[7].coords, &[0.0, ft], 1e-15);
465        array_approx_eq(&res.points[8].coords, &[et, 0.0], 1e-15);
466        array_approx_eq(&res.points[9].coords, &[ft, ft], 1e-15);
467        array_approx_eq(&res.points[10].coords, &[0.0, -4.0], 1e-15);
468        array_approx_eq(&res.points[11].coords, &[0.0, -ft], 1e-15);
469        array_approx_eq(&res.points[12].coords, &[ft, -et], 1e-15);
470        array_approx_eq(&res.points[13].coords, &[0.0, -et], 1e-15);
471        array_approx_eq(&res.points[14].coords, &[et, -ft], 1e-15);
472        array_approx_eq(&res.points[15].coords, &[ft, -ft], 1e-15);
473
474        assert_eq!(res.points[0].marker, 0);
475        assert_eq!(res.points[1].marker, 0);
476        assert_eq!(res.points[2].marker, 0);
477        assert_eq!(res.points[3].marker, -20);
478        assert_eq!(res.points[4].marker, -10);
479        assert_eq!(res.points[5].marker, 0);
480        assert_eq!(res.points[6].marker, -20);
481        assert_eq!(res.points[7].marker, -10);
482        assert_eq!(res.points[8].marker, 0);
483        assert_eq!(res.points[9].marker, 0);
484        assert_eq!(res.points[10].marker, 0);
485        assert_eq!(res.points[11].marker, -10);
486        assert_eq!(res.points[12].marker, -30);
487        assert_eq!(res.points[13].marker, -10);
488        assert_eq!(res.points[14].marker, -30);
489        assert_eq!(res.points[15].marker, 0);
490    }
491
492    #[test]
493    fn convert_tri6_to_tri3_works() {
494        #[rustfmt::skip]
495        let mesh = Mesh {
496            ndim: 2,
497            points: vec![
498                Point { id: 0, marker: 0, coords: vec![0.0, 0.0] },
499                Point { id: 1, marker: 0, coords: vec![4.0, 0.0] },
500                Point { id: 2, marker: 0, coords: vec![0.0, 4.0] },
501                Point { id: 3, marker: 0, coords: vec![2.0, 0.0] },
502                Point { id: 4, marker: -20, coords: vec![2.0, 2.0] },
503                Point { id: 5, marker: -10, coords: vec![0.0, 2.0] },
504                Point { id: 6, marker: 0, coords: vec![0.0,-4.0] },
505                Point { id: 7, marker: -30, coords: vec![2.0,-2.0] },
506                Point { id: 8, marker: -10, coords: vec![0.0,-2.0] },
507            ],
508            cells: vec![
509                Cell { id: 0, marker: 1, kind: GeoKind::Tri6, points: vec![1, 2, 0, 4, 5, 3] },
510                Cell { id: 1, marker: 2, kind: GeoKind::Tri6, points: vec![0, 6, 1, 8, 7, 3] },
511            ],
512            marked_edges: Vec::new(),
513            marked_faces: Vec::new(),
514        };
515        mesh.check_all().unwrap();
516
517        let res = mesh.convert_2d(GeoKind::Tri3).unwrap();
518
519        if SAVE_FIGURE {
520            draw(&res, false, "/tmp/gemlab/test_convert_2d_tri6_to_tri3_after.svg");
521        }
522
523        res.check_all().unwrap();
524        res.check_overlapping_points(0.01).unwrap();
525
526        assert_eq!(res.points.len(), 4);
527        assert_eq!(res.cells[0].points, &[0, 1, 2]);
528        assert_eq!(res.cells[1].points, &[2, 3, 0]);
529    }
530
531    #[test]
532    fn convert_tri3_to_tri6_works() {
533        let mesh = Samples::two_tri3().clone();
534        let res = mesh.convert_2d(GeoKind::Tri6).unwrap();
535        if SAVE_FIGURE {
536            draw(&res, false, "/tmp/gemlab/test_convert_2d_tri3_to_tri6_after.svg");
537        }
538        res.check_all().unwrap();
539        res.check_overlapping_points(0.001).unwrap();
540        assert_eq!(res.points.len(), 9);
541        assert_eq!(res.cells[0].points, (0..6).collect::<Vec<_>>());
542        assert_eq!(res.cells[1].points, &[6, 2, 1, 7, 4, 8]);
543    }
544
545    #[test]
546    fn convert_four_tri3_to_tri6_works() {
547        let mesh = Samples::four_tri3().clone();
548        let res = mesh.convert_2d(GeoKind::Tri6).unwrap();
549        if SAVE_FIGURE {
550            draw(&res, false, "/tmp/gemlab/test_convert_2d_four_tri3_to_tri6_after.svg");
551        }
552        res.check_all().unwrap();
553        res.check_overlapping_points(0.001).unwrap();
554        assert_eq!(res.points.len(), 13);
555        assert_eq!(res.cells[0].points, (0..6).collect::<Vec<_>>());
556        assert_eq!(res.cells[1].points, &[0, 6, 7, 8, 9, 10]);
557        assert_eq!(res.cells[2].points, &[0, 2, 6, 5, 11, 8]);
558        assert_eq!(res.cells[3].points, &[0, 7, 1, 10, 12, 3]);
559    }
560
561    #[test]
562    fn convert_tri3_to_tri10_works() {
563        let mesh = Samples::two_tri3().clone();
564        let res = mesh.convert_2d(GeoKind::Tri10).unwrap();
565        if SAVE_FIGURE {
566            draw(&res, false, "/tmp/gemlab/test_convert_2d_tri3_to_tri10_after.svg");
567        }
568        res.check_all().unwrap();
569        res.check_overlapping_points(0.001).unwrap();
570        assert_eq!(res.points.len(), 16);
571        assert_eq!(res.cells[0].points, (0..10).collect::<Vec<_>>());
572        assert_eq!(res.cells[1].points, &[10, 2, 1, 11, 7, 12, 13, 4, 14, 15]);
573    }
574
575    #[test]
576    fn convert_tri6_arrow_to_tri10_works() {
577        let mesh = Samples::three_tri6_arrow().clone();
578        let res = mesh.convert_2d(GeoKind::Tri10).unwrap();
579        if SAVE_FIGURE {
580            draw(&res, false, "/tmp/gemlab/test_convert_2d_tri6_arrow_to_tri10_after.svg");
581        }
582        res.check_all().unwrap();
583        res.check_overlapping_points(0.001).unwrap();
584        assert_eq!(res.points.len(), 19);
585        assert_eq!(res.cells[0].points, (0..10).collect::<Vec<_>>());
586        assert_eq!(res.cells[1].points, &[2, 1, 10, 7, 11, 12, 4, 13, 14, 15]);
587        assert_eq!(res.cells[2].points, &[0, 2, 10, 8, 14, 16, 5, 12, 17, 18]);
588        assert_eq!(res.points[0].marker, -1);
589        assert_eq!(res.points[1].marker, -2);
590        assert_eq!(res.points[2].marker, -3);
591        assert_eq!(res.points[3].marker, -10);
592        assert_eq!(res.points[4].marker, -5);
593        assert_eq!(res.points[5].marker, -6);
594        assert_eq!(res.points[6].marker, -10);
595        assert_eq!(res.points[7].marker, -5);
596        assert_eq!(res.points[8].marker, -6);
597        assert_eq!(res.points[9].marker, 0);
598        assert_eq!(res.points[10].marker, -7);
599        assert_eq!(res.points[11].marker, -20);
600        assert_eq!(res.points[12].marker, 0);
601        assert_eq!(res.points[13].marker, -20);
602        assert_eq!(res.points[14].marker, 0);
603        assert_eq!(res.points[15].marker, 0);
604        assert_eq!(res.points[16].marker, -30);
605        assert_eq!(res.points[17].marker, -30);
606        assert_eq!(res.points[18].marker, 0);
607    }
608
609    #[test]
610    fn convert_qua4_to_qua8_works() {
611        let mesh = Samples::two_qua4().clone();
612        let res = mesh.convert_2d(GeoKind::Qua8).unwrap();
613        if SAVE_FIGURE {
614            draw(&res, false, "/tmp/gemlab/test_convert_2d_qua4_to_qua8_after.svg");
615        }
616        res.check_all().unwrap();
617        res.check_overlapping_points(0.002).unwrap();
618        assert_eq!(res.points.len(), 13);
619        assert_eq!(res.cells[0].points, (0..8).collect::<Vec<_>>());
620        assert_eq!(res.cells[1].points, &[1, 8, 9, 2, 10, 11, 12, 5]);
621    }
622
623    #[test]
624    fn convert_qua12_to_qua16_works() {
625        let mesh = Samples::block_2d_four_qua12().clone();
626        let res = mesh.convert_2d(GeoKind::Qua16).unwrap();
627        if SAVE_FIGURE {
628            draw(&res, false, "/tmp/gemlab/test_convert_2d_qua12_to_qua16_after.svg");
629        }
630        res.check_all().unwrap();
631        res.check_overlapping_points(0.002).unwrap();
632        assert_eq!(res.points.len(), 33 + 4 * 4);
633        assert_eq!(res.cells[0].points, (0..16).collect::<Vec<_>>());
634        assert_eq!(
635            res.cells[1].points,
636            [1, 16, 17, 2, 18, 19, 20, 9, 21, 22, 23, 5, 24, 25, 26, 27]
637        );
638        assert_eq!(
639            res.cells[2].points,
640            [3, 2, 28, 29, 10, 30, 31, 32, 6, 33, 34, 35, 36, 37, 38, 39]
641        );
642        assert_eq!(
643            res.cells[3].points,
644            [2, 17, 40, 28, 23, 41, 42, 33, 20, 43, 44, 30, 45, 46, 47, 48]
645        );
646    }
647
648    #[test]
649    fn convert_qua12_to_qua17_works() {
650        let mesh = Samples::block_2d_four_qua12().clone();
651        let res = mesh.convert_2d(GeoKind::Qua17).unwrap();
652        if SAVE_FIGURE {
653            draw(&res, false, "/tmp/gemlab/test_convert_2d_qua12_to_qua17_after.svg");
654        }
655        res.check_all().unwrap();
656        res.check_overlapping_points(0.002).unwrap();
657        assert_eq!(res.points.len(), 3 * 9 + 6 * 3 + 4);
658        assert_eq!(res.cells[0].points, (0..17).collect::<Vec<_>>());
659        assert_eq!(
660            res.cells[1].points,
661            [1, 17, 18, 2, 19, 20, 21, 5, 22, 23, 24, 25, 26, 27, 28, 12, 11]
662        );
663        assert_eq!(
664            res.cells[2].points,
665            [3, 2, 29, 30, 6, 31, 32, 33, 34, 14, 13, 35, 36, 37, 38, 39, 40]
666        );
667        assert_eq!(
668            res.cells[3].points,
669            [2, 18, 41, 29, 21, 42, 43, 31, 44, 28, 27, 45, 46, 47, 48, 36, 35]
670        );
671        let correct = Samples::block_2d_four_qua17();
672        for i in 0..correct.cells.len() {
673            assert_eq!(&res.cells[i].points, &correct.cells[i].points);
674        }
675        let sf = 3.0 / 4.0;
676        for i in 0..correct.points.len() {
677            let scaled = &[sf * correct.points[i].coords[0], sf * correct.points[i].coords[1]];
678            array_approx_eq(&res.points[i].coords, scaled, 1e-15);
679        }
680    }
681
682    #[test]
683    fn convert_qua17_to_qua4_works() {
684        let mesh = Samples::block_2d_four_qua17().clone();
685        let res = mesh.convert_2d(GeoKind::Qua4).unwrap();
686        if SAVE_FIGURE {
687            draw(&res, false, "/tmp/gemlab/test_convert_2d_qua17_to_qua4_after.svg");
688        }
689        res.check_all().unwrap();
690        res.check_overlapping_points(0.002).unwrap();
691        assert_eq!(res.points.len(), 9);
692        let correct = Samples::block_2d_four_qua4();
693        for i in 0..correct.cells.len() {
694            assert_eq!(&res.cells[i].points, &correct.cells[i].points);
695        }
696        let sf = 4.0 / 2.0;
697        for i in 0..correct.points.len() {
698            let scaled = &[sf * correct.points[i].coords[0], sf * correct.points[i].coords[1]];
699            array_approx_eq(&res.points[i].coords, scaled, 1e-15);
700        }
701    }
702}