del_msh_cpu/
elem2elem.rs

1//! methods that generate the elements adjacent to an element
2use num_traits::AsPrimitive;
3
4pub fn face2node_of_polygon_element(num_node: usize) -> (Vec<usize>, Vec<usize>) {
5    let mut face2idx = vec![0; num_node + 1];
6    let mut idx2node = vec![0; num_node * 2];
7    for i_edge in 0..num_node {
8        face2idx[i_edge + 1] = (i_edge + 1) * 2;
9        idx2node[i_edge * 2] = i_edge;
10        idx2node[i_edge * 2 + 1] = (i_edge + 1) % num_node;
11    }
12    (face2idx, idx2node)
13}
14
15pub fn face2node_of_simplex_element(num_node: usize) -> (Vec<usize>, Vec<usize>) {
16    match num_node {
17        2 => {
18            use del_geo_core::edge::{EDGE_FACE2IDX, EDGE_IDX2NODE};
19            (EDGE_FACE2IDX.to_vec(), EDGE_IDX2NODE.to_vec())
20        }
21        3 => {
22            use del_geo_core::tri::{FACE2IDX, IDX2NODE};
23            (FACE2IDX.to_vec(), IDX2NODE.to_vec())
24        }
25        4 => (
26            del_geo_core::tet::FACE2IDX.to_vec(),
27            del_geo_core::tet::IDX2NODE.to_vec(),
28        ),
29        _ => {
30            panic!()
31        }
32    }
33}
34
35/// element adjacency of uniform mesh
36/// * `elem2vtx` - vertex index of elements
37/// * `num_node` - number of nodes par element
38/// * `vtx2elem_idx` - jagged array index of element surrounding point
39/// * `vtx2elem` - jagged array value of  element surrounding point
40///
41///  triangle: `face2jdx` = \[0,2,4,6]; `jdx2node` = \[1,2,2,0,0,1];
42pub fn from_uniform_mesh_with_vtx2elem<Index>(
43    elem2vtx: &[Index],
44    num_node: usize,
45    vtx2idx: &[Index],
46    idx2elem: &[Index],
47    face2jdx: &[usize],
48    jdx2node: &[usize],
49) -> Vec<Index>
50where
51    Index: num_traits::PrimInt + num_traits::AsPrimitive<usize>,
52    usize: num_traits::AsPrimitive<Index>,
53{
54    assert!(!vtx2idx.is_empty());
55    let num_vtx = vtx2idx.len() - 1;
56    let num_face_par_elem = face2jdx.len() - 1;
57    let num_max_node_on_face = {
58        let mut n0 = 0_usize;
59        for i_face in 0..num_face_par_elem {
60            let nno = face2jdx[i_face + 1] - face2jdx[i_face];
61            n0 = if nno > n0 { nno } else { n0 }
62        }
63        n0
64    };
65
66    let num_elem = elem2vtx.len() / num_node;
67    let mut elem2elem = vec![Index::max_value(); num_elem * num_face_par_elem];
68
69    let mut vtx2flag = vec![0; num_vtx]; // vertex index -> flag
70    let mut jdx2vtx = vec![0; num_max_node_on_face]; // face node index -> vertex index
71    for i_elem in 0..num_elem {
72        for i_face in 0..num_face_par_elem {
73            for jdx0 in 0..face2jdx[i_face + 1] - face2jdx[i_face] {
74                let i_node0 = jdx2node[jdx0 + face2jdx[i_face]];
75                assert!(i_node0 < num_node);
76                let i_vtx: usize = elem2vtx[i_elem * num_node + i_node0].as_();
77                assert!(i_vtx < num_vtx);
78                jdx2vtx[jdx0] = i_vtx;
79                vtx2flag[i_vtx] = 1;
80            }
81            let i_vtx0 = jdx2vtx[0];
82            let mut flag0 = false;
83            for &j_elem0 in &idx2elem[vtx2idx[i_vtx0].as_()..vtx2idx[i_vtx0 + 1].as_()] {
84                let j_elem0: usize = j_elem0.as_();
85                if j_elem0 == i_elem {
86                    continue;
87                }
88                for j_face in 0..num_face_par_elem {
89                    flag0 = true;
90                    for &j_node0 in &jdx2node[face2jdx[j_face]..face2jdx[j_face + 1]] {
91                        let j_vtx0: usize = elem2vtx[j_elem0 * num_node + j_node0].as_();
92                        if vtx2flag[j_vtx0] == 0 {
93                            flag0 = false;
94                            break;
95                        }
96                    }
97                    if flag0 {
98                        elem2elem[i_elem * num_face_par_elem + i_face] = j_elem0.as_();
99                        break;
100                    }
101                }
102                if flag0 {
103                    break;
104                }
105            }
106            if !flag0 {
107                elem2elem[i_elem * num_face_par_elem + i_face] = Index::max_value();
108            }
109            for ifano in 0..face2jdx[i_face + 1] - face2jdx[i_face] {
110                vtx2flag[jdx2vtx[ifano]] = 0;
111            }
112        }
113    }
114    elem2elem
115}
116
117/// element surrounding element
118/// * `elem2vtx` - vertex index of elements
119/// * `num_node` - number of nodes par element
120/// * `num_vtx` - number of vertices
121///
122///  triangle: face2idx = \[0,2,4,6]; idx2node = \[1,2,2,0,0,1];
123pub fn from_uniform_mesh<Index>(
124    elem2vtx: &[Index],
125    num_node: usize,
126    face2idx: &[usize],
127    idx2node: &[usize],
128    num_vtx: usize,
129) -> Vec<Index>
130where
131    Index: num_traits::PrimInt + num_traits::AsPrimitive<usize>,
132    usize: num_traits::AsPrimitive<Index>,
133{
134    let vtx2elem = crate::vtx2elem::from_uniform_mesh(elem2vtx, num_node, num_vtx);
135    from_uniform_mesh_with_vtx2elem(
136        elem2vtx,
137        num_node,
138        &vtx2elem.0,
139        &vtx2elem.1,
140        face2idx,
141        idx2node,
142    )
143}
144
145pub fn from_polygon_mesh_with_vtx2elem(
146    elem2idx: &[usize],
147    idx2vtx: &[usize],
148    vtx2jdx: &[usize],
149    jdx2elem: &[usize],
150) -> Vec<usize> {
151    assert!(!vtx2jdx.is_empty());
152    let num_elem = elem2idx.len() - 1;
153    let mut idx2elem = vec![usize::MAX; idx2vtx.len()];
154    for i_elem in 0..num_elem {
155        let num_edge_in_i_elem = elem2idx[i_elem + 1] - elem2idx[i_elem];
156        for i_edge in 0..num_edge_in_i_elem {
157            let i_edge2vtx = [
158                idx2vtx[elem2idx[i_elem] + i_edge],
159                idx2vtx[elem2idx[i_elem] + (i_edge + 1) % num_edge_in_i_elem],
160            ];
161            let i_vtx0 = i_edge2vtx[0];
162            for &j_elem0 in &jdx2elem[vtx2jdx[i_vtx0]..vtx2jdx[i_vtx0 + 1]] {
163                if j_elem0 == i_elem {
164                    continue;
165                }
166                let num_edge_in_j_elem0 = elem2idx[j_elem0 + 1] - elem2idx[j_elem0];
167                for j_edge in 0..num_edge_in_j_elem0 {
168                    let j_edge2vtx = [
169                        idx2vtx[elem2idx[j_elem0] + j_edge],
170                        idx2vtx[elem2idx[j_elem0] + (j_edge + 1) % num_edge_in_j_elem0],
171                    ];
172                    if i_edge2vtx[0] != j_edge2vtx[1] || i_edge2vtx[1] != j_edge2vtx[0] {
173                        continue;
174                    }
175                    idx2elem[elem2idx[i_elem] + i_edge] = j_elem0;
176                    break;
177                }
178                if idx2elem[elem2idx[i_elem] + i_edge] != usize::MAX {
179                    break;
180                }
181            }
182        }
183    }
184    idx2elem
185}
186
187pub fn from_polygon_mesh(elem2idx: &[usize], idx2vtx: &[usize], num_vtx: usize) -> Vec<usize> {
188    let vtx2elem = crate::vtx2elem::from_polygon_mesh(elem2idx, idx2vtx, num_vtx);
189    from_polygon_mesh_with_vtx2elem(elem2idx, idx2vtx, &vtx2elem.0, &vtx2elem.1)
190}
191
192/// Extract the boundary surface mesh from a uniform volumetric mesh.
193///
194/// A face is on the boundary when its `elem2elem` entry equals `Index::max_value()`.
195/// Returns a flat array of vertex indices for the boundary faces (uniform surface mesh).
196/// The number of nodes per boundary face is `face2idx[1] - face2idx[0]`.
197///
198/// # Arguments
199/// * `elem2vtx` - vertex indices of elements, length `num_elem * num_node`
200/// * `num_node` - number of nodes per element (e.g. 4 for tets)
201/// * `elem2elem` - element adjacency array, length `num_elem * num_face_per_elem`;
202///   boundary faces have value `Index::max_value()`
203/// * `face2idx` - CSR offsets into `idx2node` for each face of an element
204/// * `idx2node` - local node indices on each face
205pub fn extract_boundary_mesh_for_uniform_mesh<Index>(
206    elem2vtx: &[Index],
207    num_node: usize,
208    elem2elem: &[Index],
209    face2idx: &[usize],
210    idx2node: &[usize],
211) -> Vec<Index>
212where
213    Index: num_traits::PrimInt + num_traits::AsPrimitive<usize>,
214    usize: num_traits::AsPrimitive<Index>,
215{
216    let num_face_per_elem = face2idx.len() - 1;
217    let num_elem = elem2vtx.len() / num_node;
218    let mut bnd_face2vtx = Vec::<Index>::new();
219    for i_elem in 0..num_elem {
220        for i_face in 0..num_face_per_elem {
221            if elem2elem[i_elem * num_face_per_elem + i_face] != Index::max_value() {
222                continue;
223            }
224            #[allow(clippy::needless_range_loop)]
225            for jdx in face2idx[i_face]..face2idx[i_face + 1] {
226                let i_node = idx2node[jdx];
227                bnd_face2vtx.push(elem2vtx[i_elem * num_node + i_node]);
228            }
229        }
230    }
231    bnd_face2vtx
232}
del_msh_cpu/elem2elem.rs

del_msh_cpu/
elem2elem.rs
