del_msh_cpu/
polyloop.rs

1//! methods for 2D and 3D poly loop
2
3use num_traits::AsPrimitive;
4
5/// return  arc-length of a 2D or 3D poly loop
6pub fn arclength_from_vtx2vecn<T, const N: usize>(vtxs: &[[T; N]]) -> T
7where
8    T: num_traits::Float + Copy + 'static + std::iter::Sum,
9    f64: AsPrimitive<T>,
10{
11    use del_geo_core::vecn::VecN;
12    if vtxs.len() < 2 {
13        return T::zero();
14    }
15    let np = vtxs.len();
16    let mut len: T = T::zero();
17    for ip0 in 0..np {
18        let ip1 = (ip0 + 1) % np;
19        len = len + vtxs[ip0].sub(&vtxs[ip1]).norm();
20    }
21    len
22}
23
24/// return  arc-length of a 2D or 3D poly loop
25pub fn arclength<T, const N: usize>(vtx2xyz: &[T]) -> T
26where
27    T: num_traits::Float,
28{
29    let np = vtx2xyz.len() / N;
30    let mut len: T = T::zero();
31    for ip0 in 0..np {
32        let ip1 = (ip0 + 1) % np;
33        let p0: &[T; N] = &vtx2xyz[ip0 * N..ip0 * N + N].try_into().unwrap();
34        let p1: &[T; N] = &vtx2xyz[ip1 * N..ip1 * N + N].try_into().unwrap();
35        len = len + del_geo_core::edge::length::<T, N>(p0, p1);
36    }
37    len
38}
39
40pub fn edge2length<T, const N: usize>(vtx2xyz: &[T]) -> Vec<T>
41where
42    T: num_traits::Float + std::ops::AddAssign,
43{
44    let np = vtx2xyz.len() / N;
45    let mut edge2length = Vec::<T>::with_capacity(np);
46    for ip0 in 0..np {
47        let ip1 = (ip0 + 1) % np;
48        let p0: &[T; N] = &vtx2xyz[ip0 * N..ip0 * N + N].try_into().unwrap();
49        let p1: &[T; N] = &vtx2xyz[ip1 * N..ip1 * N + N].try_into().unwrap();
50        edge2length.push(del_geo_core::edge::length::<T, N>(p0, p1));
51    }
52    edge2length
53}
54
55/// the center of gravity for polyloop.
56/// Here polyloop is a looped wire, not the polygonal face bounded by the polyloop
57pub fn cog_as_edges<T, const N: usize>(vtx2xyz: &[T]) -> [T; N]
58where
59    T: num_traits::Float + Copy + 'static + std::iter::Sum,
60    f64: AsPrimitive<T>,
61{
62    let num_vtx = vtx2xyz.len() / N;
63    assert_eq!(vtx2xyz.len(), num_vtx * N);
64    let mut cog = [T::zero(); N];
65    let mut len = T::zero();
66    use del_geo_core::vecn::VecN;
67    for i_edge in 0..num_vtx {
68        let iv0 = i_edge;
69        let iv1 = (i_edge + 1) % num_vtx;
70        let q0: &[T; N] = vtx2xyz[iv0 * N..iv0 * N + N].try_into().unwrap();
71        let q1: &[T; N] = vtx2xyz[iv1 * N..iv1 * N + N].try_into().unwrap();
72        let l = q0.sub(q1).norm();
73        let d = q0.add(q1).scale(0.5_f64.as_() * l);
74        cog = cog.add(&d);
75        len = len + l;
76    }
77    cog.scale(T::one() / len)
78}
79
80#[test]
81fn test_cog() {
82    let mut vtx2xy = crate::polyloop2::from_circle(1f32, 32);
83    let (x0, y0) = (1.3, 0.5);
84    vtx2xy.chunks_mut(2).for_each(|v| {
85        v[0] += x0;
86        v[1] += y0
87    });
88    let cog = cog_as_edges::<f32, 2>(vtx2xy.as_slice());
89    assert!(
90        del_geo_core::edge::length::<f32, 2>(&[x0, y0], cog.as_slice().try_into().unwrap())
91            < 1.0e-5
92    );
93}
94
95/*
96/// the center of gravity for polyloop.
97/// Here "polyloop" is a looped wire, not the polygonal face bounded by the polyloop
98pub fn cog_from_vtx2vecn_as_edges<T, const N: usize>(
99    vtx2xyz: &[nalgebra::SVector<T, N>],
100) -> [T; N]
101where
102    T: num_traits::Float + Copy + 'static,
103    f64: AsPrimitive<T>,
104{
105    use del_geo_core::vecn::Arr;
106    let num_vtx = vtx2xyz.len() / 3;
107    assert_eq!(vtx2xyz.len(), num_vtx);
108    let mut cog = [T::zero(); N];
109    let mut len = T::zero();
110    for i_edge in 0..num_vtx {
111        let iv0 = i_edge;
112        let iv1 = (i_edge + 1) % num_vtx;
113        let q0 = &vtx2xyz[iv0];
114        let q1 = &vtx2xyz[iv1];
115        let l = q0.sub(q1).norm();
116        let d = (q0 + q1).scale(0.5_f64.as_() * l);
117        cog = cog.add(&d);
118        len += l;
119    }
120    cog.scale(T::one()/len)
121}
122 */
123
124pub fn cov<T, const N: usize>(vtx2xyz: &[T]) -> [[T; N]; N]
125where
126    T: num_traits::Float + Copy + 'static + std::iter::Sum,
127    f64: AsPrimitive<T>,
128{
129    use del_geo_core::vecn::VecN;
130    let one = T::one();
131    let three = one + one + one;
132    let six = three + three;
133    let num_vtx = vtx2xyz.len() / N;
134    assert_eq!(vtx2xyz.len(), num_vtx * N);
135    let cog = cog_as_edges::<T, N>(vtx2xyz);
136    let mut cov = [[T::zero(); N]; N];
137    for i_edge in 0..num_vtx {
138        let iv0 = i_edge;
139        let iv1 = (i_edge + 1) % num_vtx;
140        let q0 = crate::vtx2xn::to_xn(vtx2xyz, iv0).sub(&cog);
141        let q1 = crate::vtx2xn::to_xn(vtx2xyz, iv1).sub(&cog);
142        let l = q0.sub(&q1).norm();
143        for i in 0..N {
144            for j in 0..N {
145                cov[i][j] = cov[i][j]
146                    + (q0[i] * q0[j] + q1[i] * q1[j]) * (l / three)
147                    + (q0[i] * q1[j] + q1[i] * q0[j]) * (l / six);
148            }
149        }
150        // cov += (q0 * q0.transpose() + q1 * q1.transpose()).scale(l / 3_f64.as_());
151        //cov += (q0 * q1.transpose() + q1 * q0.transpose()).scale(l / 6_f64.as_());
152    }
153    cov
154}
155
156pub fn resample<T, const N: usize>(vtx2xyz_in: &[T], num_edge_out: usize) -> Vec<T>
157where
158    T: num_traits::Float + Copy + std::iter::Sum + 'static,
159    usize: AsPrimitive<T>,
160{
161    use del_geo_core::vecn::VecN;
162    let mut v2x_out = Vec::<T>::new();
163    let num_edge_in = vtx2xyz_in.len() / N;
164    let len_edge_out = arclength::<T, N>(vtx2xyz_in) / num_edge_out.as_();
165    v2x_out.extend_from_slice(&vtx2xyz_in[0..N]);
166    let mut i_edge_in = 0;
167    let mut traveled_ratio0 = T::zero();
168    let mut remaining_length = len_edge_out;
169    loop {
170        if i_edge_in >= num_edge_in {
171            break;
172        }
173        if v2x_out.len() >= num_edge_out * N {
174            break;
175        }
176        let i0 = i_edge_in;
177        let i1 = (i_edge_in + 1) % num_edge_in;
178        let p0: &[T; N] = vtx2xyz_in[i0 * N..i0 * N + N].try_into().unwrap();
179        let p1: &[T; N] = vtx2xyz_in[i1 * N..i1 * N + N].try_into().unwrap();
180        let len_edge0 = p1.sub(p0).norm();
181        let len_togo0 = len_edge0 * (T::one() - traveled_ratio0);
182        if len_togo0 > remaining_length {
183            // put point in this segment
184            traveled_ratio0 = traveled_ratio0 + remaining_length / len_edge0;
185            let pn = p0
186                .scale(T::one() - traveled_ratio0)
187                .add(&p1.scale(traveled_ratio0));
188            v2x_out.extend(pn.iter());
189            remaining_length = len_edge_out;
190        } else {
191            // next segment
192            remaining_length = remaining_length - len_togo0;
193            traveled_ratio0 = T::zero();
194            i_edge_in += 1;
195        }
196    }
197    v2x_out
198}
199
200pub fn resample_multiple_loops_remain_original_vtxs<T>(
201    loop2idx_inout: &mut Vec<usize>,
202    idx2vtx_inout: &mut Vec<usize>,
203    vtx2vec_inout: &mut Vec<[T; 2]>,
204    max_edge_length: T,
205) where
206    T: num_traits::Float + Copy + AsPrimitive<usize>,
207    usize: AsPrimitive<T>,
208{
209    use del_geo_core::vec2::Vec2;
210    let one = T::one();
211    assert_eq!(vtx2vec_inout.len(), idx2vtx_inout.len());
212    let loop2idx_in = loop2idx_inout.clone();
213    let idx2vtx_in = idx2vtx_inout.clone();
214    assert!(idx2vtx_in.len() >= 2);
215    let num_loop = loop2idx_in.len() - 1;
216    let mut edge2point: Vec<Vec<usize>> = vec![vec!(); idx2vtx_in.len()];
217    {
218        for i_loop in 0..num_loop {
219            assert!(loop2idx_in[i_loop + 1] > loop2idx_in[i_loop]);
220            let np = loop2idx_in[i_loop + 1] - loop2idx_in[i_loop];
221            for ip in 0..np {
222                let iipo0 = loop2idx_in[i_loop] + ip;
223                let iipo1 = loop2idx_in[i_loop] + (ip + 1) % np;
224                assert!(iipo0 < idx2vtx_in.len());
225                assert!(iipo1 < idx2vtx_in.len());
226                let ipo0 = idx2vtx_in[iipo0];
227                let ipo1 = idx2vtx_in[iipo1];
228                assert!(ipo0 < vtx2vec_inout.len());
229                assert!(ipo1 < vtx2vec_inout.len());
230                let po0 = vtx2vec_inout[ipo0]; // never use reference here because aVec2 will resize afterward
231                let po1 = vtx2vec_inout[ipo1]; // never use reference here because aVec2 will resize afterward
232                let nadd: usize = (po0.sub(&po1).norm() / max_edge_length).as_();
233                if nadd == 0 {
234                    continue;
235                }
236                for iadd in 0..nadd {
237                    let r2: T = (iadd + 1).as_() / (nadd + 1).as_();
238                    let v2 = po0.scale(one - r2).add(&po1.scale(r2));
239                    let ipo2 = vtx2vec_inout.len();
240                    vtx2vec_inout.push(v2);
241                    assert!(iipo0 < edge2point.len());
242                    edge2point[iipo0].push(ipo2);
243                }
244            }
245        }
246    }
247    ////
248    loop2idx_inout.resize(num_loop + 1, usize::MAX);
249    loop2idx_inout[0] = 0;
250    for iloop in 0..num_loop {
251        let nbar0 = loop2idx_in[iloop + 1] - loop2idx_in[iloop];
252        let mut nbar1 = nbar0;
253        for ibar in 0..nbar0 {
254            let iip_loop = loop2idx_in[iloop] + ibar;
255            nbar1 += edge2point[iip_loop].len();
256        }
257        loop2idx_inout[iloop + 1] = loop2idx_inout[iloop] + nbar1;
258    }
259    // adding new vertices on the outline
260    idx2vtx_inout.resize(loop2idx_inout[num_loop], usize::MAX);
261    let mut i_vtx0 = 0;
262    for i_loop in 0..num_loop {
263        for iip_loop in loop2idx_in[i_loop]..loop2idx_in[i_loop + 1] {
264            let ip_loop = idx2vtx_in[iip_loop];
265            idx2vtx_inout[i_vtx0] = ip_loop;
266            i_vtx0 += 1;
267            #[allow(clippy::needless_range_loop)]
268            for iadd in 0..edge2point[ip_loop].len() {
269                idx2vtx_inout[i_vtx0] = edge2point[iip_loop][iadd];
270                i_vtx0 += 1;
271            }
272        }
273    }
274    assert_eq!(idx2vtx_inout.len(), vtx2vec_inout.len());
275    assert_eq!(idx2vtx_inout.len(), i_vtx0);
276}
277
278pub fn to_cylinder_trimeshes<Real>(
279    vtx2xy: &[Real],
280    num_dim: usize,
281    radius: Real,
282) -> (Vec<usize>, Vec<Real>)
283where
284    Real: num_traits::Float + num_traits::FloatConst + 'static + Copy,
285    usize: AsPrimitive<Real>,
286{
287    let num_vtx = vtx2xy.len() / num_dim;
288    let mut out_tri2vtx: Vec<usize> = vec![];
289    let mut out_vtx2xyz: Vec<Real> = vec![];
290    for i_edge in 0..num_vtx {
291        let i0 = i_edge;
292        let i1 = (i_edge + 1) % num_vtx;
293        let p0 = &vtx2xy[i0 * num_dim..(i0 + 1) * num_dim];
294        let p1 = &vtx2xy[i1 * num_dim..(i1 + 1) * num_dim];
295        let p0 = if num_dim == 3 {
296            [p0[0], p0[1], p0[2]]
297        } else {
298            [p0[0], p0[1], Real::zero()]
299        };
300        let p1 = if num_dim == 3 {
301            [p1[0], p1[1], p1[2]]
302        } else {
303            [p1[0], p1[1], Real::zero()]
304        };
305        let (tri2vtx, vtx2xyz) =
306            crate::trimesh3_primitive::cylinder_open_connecting_two_points(32, radius, &p0, &p1);
307        crate::uniform_mesh::merge(
308            &mut out_tri2vtx,
309            &mut out_vtx2xyz,
310            tri2vtx.as_slice(),
311            vtx2xyz.as_slice(),
312            3,
313        );
314    }
315    (out_tri2vtx, out_vtx2xyz)
316}
317
318pub fn edge2vtx(num_vtx: usize) -> Vec<usize> {
319    let mut edge2vtx = Vec::<usize>::with_capacity(num_vtx * 2);
320    for i_vtx in 0..num_vtx {
321        edge2vtx.push(i_vtx);
322        edge2vtx.push((i_vtx + 1) % num_vtx);
323    }
324    edge2vtx
325}
326
327pub fn flip<Real>(vtx2xy: &[Real], num_dim: usize) -> Vec<Real>
328where
329    Real: num_traits::Float,
330{
331    let mut vtx2xy_out = Vec::<Real>::with_capacity(vtx2xy.len());
332    for xy in vtx2xy.chunks(num_dim).rev() {
333        vtx2xy_out.extend(xy);
334        /*
335        for x in xy.iter() {
336            vtx2xy_out.push(*x);
337        }
338         */
339    }
340    vtx2xy_out
341}
del_msh_cpu/polyloop.rs

del_msh_cpu/
polyloop.rs
