del-msh-cpu 0.1.49

mesh utility library for computer graphics research and prototyping
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214

/// Given face `i_face` of tet `vtxs_i`, find the local face index on `j_tet` that shares the same
/// three vertices. Matching is done by vertex-index sum, so vertex order does not matter.
pub fn find_adjacent_face_index(
    vtxs_i: &[usize; 4],
    i_face: usize,
    j_tet: usize,
    tet2vtx: &[[usize; 4]],
    face2idx: &[usize],
    idx2node: &[usize],
) -> usize {
    let iv0 = vtxs_i[idx2node[face2idx[i_face]]];
    let iv1 = vtxs_i[idx2node[face2idx[i_face] + 1]];
    let iv2 = vtxs_i[idx2node[face2idx[i_face] + 2]];
    let i_sum = iv0 + iv1 + iv2;
    assert_ne!(iv0, iv1);
    assert_ne!(iv0, iv2);
    let vtxs_j = &tet2vtx[j_tet];
    for j_face in 0..4 {
        let jv0 = vtxs_j[idx2node[face2idx[j_face]]];
        let jv1 = vtxs_j[idx2node[face2idx[j_face] + 1]];
        let jv2 = vtxs_j[idx2node[face2idx[j_face] + 2]];
        let j_sum = jv0 + jv1 + jv2;
        if i_sum == j_sum {
            return j_face;
        }
    }
    dbg!(iv0, iv1, iv2, &vtxs_j);
    panic!();
}

/// Logically delete `i_tet` from the mesh: severs all adjacency links to its neighbours in
/// `tet2tet` (both sides) and overwrites its vertex slots in `tet2vtx` with `usize::MAX`.
pub fn remove(i_tet: usize, tet2tet: &mut [usize], tet2vtx: &mut [[usize; 4]]) {
    for i2_node in 0..4 {
        let k_tet = tet2tet[i_tet * 4 + i2_node];
        if k_tet == usize::MAX {
            continue;
        }
        let k2_node = find_adjacent_face_index(
            &tet2vtx[i_tet],
            i2_node,
            k_tet,
            tet2vtx,
            &del_geo_core::tet::FACE2IDX,
            &del_geo_core::tet::IDX2NODE,
        );
        tet2tet[i_tet * 4 + i2_node] = usize::MAX;
        tet2tet[k_tet * 4 + k2_node] = usize::MAX;
    }
    tet2vtx[i_tet].iter_mut().for_each(|v| *v = usize::MAX);
}

/// Candidate for collapsing two adjacent boundary tets into a single pyramid.
/// `i0_node` is the boundary face node of `i_tet`; `i1_node` / `j0_node` identify the shared
/// interior apex vertex.
#[derive(Debug, Clone)]
struct MergeTwoTetsIntoPrm {
    i_tet: usize,
    i0_node: usize,
    i1_node: usize,
    j_tet: usize,
    j0_node: usize,
}

impl MergeTwoTetsIntoPrm {
    /// Returns the 5 vertex indices of the pyramid formed by merging the two tets.
    /// Winding: [base quad (k0,k3,k2,k1), apex k4].
    fn prism_vtxs(&self, tet2vtx: &[[usize; 4]]) -> [usize; 5] {
        use del_geo_core::tet::{FACE2IDX, IDX2NODE};
        let (i_tet, j_tet, i0_node, i1_node, j0_node) = (
            self.i_tet,
            self.j_tet,
            self.i0_node,
            self.i1_node,
            self.j0_node,
        );
        let k4_vtx = tet2vtx[i_tet][i0_node];
        assert_eq!(tet2vtx[j_tet][j0_node], k4_vtx);
        let k0_vtx = tet2vtx[i_tet][i1_node];
        let i0_nofa = (0..3)
            .find(|i| IDX2NODE[FACE2IDX[i0_node] + i] == i1_node)
            .unwrap();
        let i2_node = IDX2NODE[FACE2IDX[i0_node] + (i0_nofa + 1) % 3];
        let i3_node = IDX2NODE[FACE2IDX[i0_node] + (i0_nofa + 2) % 3];
        let j1_node = find_adjacent_face_index(
            &tet2vtx[i_tet],
            i1_node,
            j_tet,
            tet2vtx,
            &FACE2IDX,
            &IDX2NODE,
        );
        let k1_vtx = tet2vtx[i_tet][i2_node];
        let k2_vtx = tet2vtx[j_tet][j1_node];
        let k3_vtx = tet2vtx[i_tet][i3_node];
        [k0_vtx, k3_vtx, k2_vtx, k1_vtx, k4_vtx]
    }
}

/// Scan boundary tet pairs and greedily merge those whose resulting pyramid has a better
/// condition number (>1.3×) than both source tets. Returns `(pyrmd2vtx, consumed_tets)` where
/// `pyrmd2vtx` is a flat list of 5-vertex pyramids and `consumed_tets` is the set of tet indices
/// that were absorbed. Each tet is used at most once.
pub fn make_pyramids_on_boundary(
    tet2vtx: &[[usize; 4]],
    tet2tet: &[usize],
    vtx2xyz: &[[f64; 3]],
) -> (Vec<[usize; 5]>, Vec<[usize; 2]>) {
    let num_tet = tet2vtx.len();
    let mut cands: Vec<MergeTwoTetsIntoPrm> = vec![];
    for i_tet in 0..num_tet {
        let cands_i: Vec<MergeTwoTetsIntoPrm> = (0..4)
            .flat_map(|i0_node| {
                if tet2tet[i_tet * 4 + i0_node] != usize::MAX {
                    return None;
                }
                let i0_vtx = tet2vtx[i_tet][i0_node];
                let res = (0..3)
                    .filter_map(|i| {
                        let i1_node = (i0_node + i) % 4;
                        let j_tet = tet2tet[i_tet * 4 + i1_node];
                        if j_tet == usize::MAX {
                            return None;
                        }
                        let j0_node = (0..4).find(|&i| tet2vtx[j_tet][i] == i0_vtx).unwrap();
                        if tet2tet[j_tet * 4 + j0_node] != usize::MAX {
                            return None;
                        }
                        Some(MergeTwoTetsIntoPrm {
                            i_tet,
                            i0_node,
                            i1_node,
                            j_tet,
                            j0_node,
                        })
                    })
                    .collect::<Vec<MergeTwoTetsIntoPrm>>();
                Some(res)
            })
            .flatten()
            .collect();
        cands.extend_from_slice(&cands_i);
    }

    let cands: Vec<_> = cands
        .into_iter()
        .filter(|cand| {
            let ci = {
                let i_tet = cand.i_tet;
                let p0i = &vtx2xyz[tet2vtx[i_tet][0]];
                let p1i = &vtx2xyz[tet2vtx[i_tet][1]];
                let p2i = &vtx2xyz[tet2vtx[i_tet][2]];
                let p3i = &vtx2xyz[tet2vtx[i_tet][3]];
                1.0 / del_geo_core::tet::condition_number(p0i, p1i, p2i, p3i).unwrap()
            };
            let cj = {
                let j_tet = cand.j_tet;
                let p0j = &vtx2xyz[tet2vtx[j_tet][0]];
                let p1j = &vtx2xyz[tet2vtx[j_tet][1]];
                let p2j = &vtx2xyz[tet2vtx[j_tet][2]];
                let p3j = &vtx2xyz[tet2vtx[j_tet][3]];
                1.0 / del_geo_core::tet::condition_number(p0j, p1j, p2j, p3j).unwrap()
            };
            let c_new = {
                let pnode2vtx = cand.prism_vtxs(tet2vtx);
                let p0 = &vtx2xyz[pnode2vtx[0]];
                let p1 = &vtx2xyz[pnode2vtx[1]];
                let p2 = &vtx2xyz[pnode2vtx[2]];
                let p3 = &vtx2xyz[pnode2vtx[3]];
                let p4 = &vtx2xyz[pnode2vtx[4]];
                use del_geo_core::pyramid::jacobian_determinant_and_conditiona_number as djac_cnd;
                let (d0, c0) = djac_cnd(p0, p1, p2, p3, p4, &[0.5, 0.5, 0.5]).unwrap();
                let (d1, c1) = djac_cnd(p0, p1, p2, p3, p4, &[0.5, 0.5, 0.0]).unwrap();
                let (d2, c2) = djac_cnd(p0, p1, p2, p3, p4, &[0.0, 0.0, 0.0]).unwrap();
                let (d3, c3) = djac_cnd(p0, p1, p2, p3, p4, &[1.0, 0.0, 0.0]).unwrap();
                let (d4, c4) = djac_cnd(p0, p1, p2, p3, p4, &[1.0, 1.0, 0.0]).unwrap();
                let (d5, c5) = djac_cnd(p0, p1, p2, p3, p4, &[0.0, 1.0, 0.0]).unwrap();
                let (c0, c1, c2, c3, c4, c5) =
                    (1.0 / c0, 1.0 / c1, 1.0 / c2, 1.0 / c3, 1.0 / c4, 1.0 / c5);
                let dmin = [d0, d1, d2, d3, d4, d5]
                    .iter()
                    .copied()
                    .min_by(|x, y| x.total_cmp(y));
                dmin.filter(|&d| d > 0.0).and_then(|_| {
                    [c0, c1, c2, c3, c4, c5]
                        .iter()
                        .copied()
                        .min_by(|x, y| x.total_cmp(y))
                })
            };
            //
            let c_old = ci.min(cj);
            c_new.is_some_and(|c| c > c_old * 1.3)
        })
        .collect();
    dbg!(cands.len());

    let mut seen = std::collections::HashSet::<usize>::new();
    let mut pyrmd2vtx: Vec<[usize; 5]> = vec![];
    let mut consumed: Vec<[usize; 2]> = vec![];
    for cand in &cands {
        let i_tet = cand.i_tet;
        let j_tet = cand.j_tet;
        if seen.contains(&i_tet) || seen.contains(&j_tet) {
            continue;
        }
        let pnode2vtx = cand.prism_vtxs(tet2vtx);
        pyrmd2vtx.push(pnode2vtx);
        consumed.push([i_tet, j_tet]);
        seen.insert(i_tet);
        seen.insert(j_tet);
    }
    (pyrmd2vtx, consumed)
}