ndgrid/geometry/
geometry_map.rs

1//! Geometry map
2use crate::{traits::GeometryMap as GeometryMapTrait, types::Scalar};
3use itertools::izip;
4use ndelement::{reference_cell, traits::MappedFiniteElement, types::ReferenceCellType};
5use rlst::{Array, DynArray, MutableArrayImpl, RlstScalar, ValueArrayImpl};
6
7/// Single element geometry
8pub struct GeometryMap<'a, T, B2D, C2D> {
9    geometry_points: &'a Array<B2D, 2>,
10    entities: &'a Array<C2D, 2>,
11    tdim: usize,
12    gdim: usize,
13    table: DynArray<T, 4>,
14}
15
16/// Dot product of two vectors
17fn dot<T: RlstScalar, Array1Impl: ValueArrayImpl<T, 1>>(
18    a: &Array<Array1Impl, 1>,
19    b: &Array<Array1Impl, 1>,
20) -> T {
21    debug_assert!(a.shape()[0] == b.shape()[0]);
22    izip!(a.iter_value(), b.iter_value())
23        .map(|(i, j)| i * j)
24        .sum::<T>()
25}
26
27/// Invert a matrix and return the determinad, or compute the Moore-Penrose left pseudoinverse
28/// and return sqrt(det(A^T*A) if the matrix is not square
29/// Note: the shape argument assumes column-major ordering
30fn inverse_and_det<
31    T: RlstScalar,
32    Array2Impl: ValueArrayImpl<T, 2>,
33    Array2ImplMut: MutableArrayImpl<T, 2>,
34>(
35    mat: &Array<Array2Impl, 2>,
36    result: &mut Array<Array2ImplMut, 2>,
37) -> T {
38    let gdim = mat.shape()[0];
39    let tdim = mat.shape()[1];
40
41    debug_assert!(result.shape()[0] == tdim);
42    debug_assert!(result.shape()[1] == gdim);
43    debug_assert!(tdim <= gdim);
44
45    match tdim {
46        0 => T::one(),
47        1 => {
48            let det = dot(&mat.r().slice(1, 0), &mat.r().slice(1, 0));
49            for (r, m) in izip!(result.iter_mut(), mat.iter_value()) {
50                *r = m / det;
51            }
52            det
53        }
54        2 => {
55            let col0 = &mat.r().slice(1, 0);
56            let col1 = &mat.r().slice(1, 1);
57            let ata = [
58                dot(col0, col0),
59                dot(col0, col1),
60                dot(col1, col0),
61                dot(col1, col1),
62            ];
63            let ata_det = ata[0] * ata[3] - ata[1] * ata[2];
64            let ata_inv = [
65                ata[3] / ata_det,
66                -ata[2] / ata_det,
67                -ata[1] / ata_det,
68                ata[0] / ata_det,
69            ];
70
71            for i in 0..gdim {
72                *result.get_mut([0, i]).unwrap() = ata_inv[0] * mat.get_value([i, 0]).unwrap()
73                    + ata_inv[2] * mat.get_value([i, 1]).unwrap();
74                *result.get_mut([1, i]).unwrap() = ata_inv[1] * mat.get_value([i, 0]).unwrap()
75                    + ata_inv[3] * mat.get_value([i, 1]).unwrap();
76            }
77            ata_det
78        }
79        3 => {
80            let col0 = &mat.r().slice(1, 0);
81            let col1 = &mat.r().slice(1, 1);
82            let col2 = &mat.r().slice(1, 2);
83            let ata = [
84                dot(col0, col0),
85                dot(col0, col1),
86                dot(col0, col2),
87                dot(col1, col0),
88                dot(col1, col1),
89                dot(col1, col2),
90                dot(col2, col0),
91                dot(col2, col1),
92                dot(col2, col2),
93            ];
94            let ata_det = ata[0] * (ata[4] * ata[8] - ata[5] * ata[7])
95                + ata[1] * (ata[5] * ata[6] - ata[3] * ata[8])
96                + ata[2] * (ata[3] * ata[7] - ata[4] * ata[6]);
97            let ata_inv = [
98                (ata[4] * ata[8] - ata[5] * ata[7]) / ata_det,
99                (ata[2] * ata[7] - ata[1] * ata[8]) / ata_det,
100                (ata[1] * ata[5] - ata[2] * ata[4]) / ata_det,
101                (ata[5] * ata[6] - ata[3] * ata[8]) / ata_det,
102                (ata[0] * ata[8] - ata[2] * ata[6]) / ata_det,
103                (ata[2] * ata[3] - ata[0] * ata[5]) / ata_det,
104                (ata[3] * ata[7] - ata[4] * ata[6]) / ata_det,
105                (ata[1] * ata[6] - ata[0] * ata[7]) / ata_det,
106                (ata[0] * ata[4] - ata[1] * ata[3]) / ata_det,
107            ];
108
109            for i in 0..gdim {
110                *result.get_mut([0, i]).unwrap() = ata_inv[0] * mat.get_value([i, 0]).unwrap()
111                    + ata_inv[3] * mat.get_value([i, 1]).unwrap()
112                    + ata_inv[6] * mat.get_value([i, 2]).unwrap();
113                *result.get_mut([1, i]).unwrap() = ata_inv[1] * mat.get_value([i, 0]).unwrap()
114                    + ata_inv[4] * mat.get_value([i, 1]).unwrap()
115                    + ata_inv[7] * mat.get_value([i, 2]).unwrap();
116                *result.get_mut([2, i]).unwrap() = ata_inv[2] * mat.get_value([i, 0]).unwrap()
117                    + ata_inv[5] * mat.get_value([i, 1]).unwrap()
118                    + ata_inv[8] * mat.get_value([i, 2]).unwrap();
119            }
120            ata_det
121        }
122        _ => {
123            panic!("Unsupported dimension");
124        }
125    }
126    .sqrt()
127}
128
129fn cross<T: RlstScalar, Array2Impl: ValueArrayImpl<T, 2>, Array1ImplMut: MutableArrayImpl<T, 1>>(
130    mat: &Array<Array2Impl, 2>,
131    result: &mut Array<Array1ImplMut, 1>,
132) {
133    debug_assert!(mat.shape()[0] == mat.shape()[1] + 1);
134    match mat.shape()[1] {
135        0 => {}
136        1 => {
137            debug_assert!(result.shape()[0] == 2);
138            *result.get_mut([0]).unwrap() = mat.get_value([1, 0]).unwrap();
139            *result.get_mut([1]).unwrap() = -mat.get_value([0, 0]).unwrap();
140        }
141        2 => {
142            debug_assert!(result.shape()[0] == 3);
143            *result.get_mut([0]).unwrap() = mat.get_value([1, 0]).unwrap()
144                * mat.get_value([2, 1]).unwrap()
145                - mat.get_value([2, 0]).unwrap() * mat.get_value([1, 1]).unwrap();
146            *result.get_mut([1]).unwrap() = mat.get_value([2, 0]).unwrap()
147                * mat.get_value([0, 1]).unwrap()
148                - mat.get_value([0, 0]).unwrap() * mat.get_value([2, 1]).unwrap();
149            *result.get_mut([2]).unwrap() = mat.get_value([0, 0]).unwrap()
150                * mat.get_value([1, 1]).unwrap()
151                - mat.get_value([1, 0]).unwrap() * mat.get_value([0, 1]).unwrap();
152        }
153        _ => {
154            unimplemented!();
155        }
156    }
157}
158
159impl<'a, T: Scalar, B2D: ValueArrayImpl<T, 2>, C2D: ValueArrayImpl<usize, 2>>
160    GeometryMap<'a, T, B2D, C2D>
161{
162    /// Create new
163    pub fn new<A2D: ValueArrayImpl<T, 2>>(
164        element: &impl MappedFiniteElement<CellType = ReferenceCellType, T = T>,
165        points: &Array<A2D, 2>,
166        geometry_points: &'a Array<B2D, 2>,
167        entities: &'a Array<C2D, 2>,
168    ) -> Self {
169        let tdim = reference_cell::dim(element.cell_type());
170        debug_assert!(points.shape()[0] == tdim);
171        let gdim = geometry_points.shape()[0];
172        let npoints = points.shape()[1];
173
174        let mut table = DynArray::<T, 4>::from_shape(element.tabulate_array_shape(1, npoints));
175        element.tabulate(points, 1, &mut table);
176
177        Self {
178            geometry_points,
179            entities,
180            tdim,
181            gdim,
182            table,
183        }
184    }
185}
186
187impl<T: Scalar, B2D: ValueArrayImpl<T, 2>, C2D: ValueArrayImpl<usize, 2>> GeometryMapTrait
188    for GeometryMap<'_, T, B2D, C2D>
189{
190    type T = T;
191
192    fn entity_topology_dimension(&self) -> usize {
193        self.tdim
194    }
195    fn geometry_dimension(&self) -> usize {
196        self.gdim
197    }
198    fn point_count(&self) -> usize {
199        self.table.shape()[1]
200    }
201    fn physical_points<Array2Impl: MutableArrayImpl<T, 2>>(
202        &self,
203        entity_index: usize,
204        points: &mut Array<Array2Impl, 2>,
205    ) {
206        let npts = self.table.shape()[1];
207        debug_assert!(points.shape()[0] == self.gdim);
208        debug_assert!(points.shape()[1] == npts);
209
210        points.fill_with_value(T::default());
211        for i in 0..self.entities.shape()[0] {
212            let v = unsafe { self.entities.get_value_unchecked([i, entity_index]) };
213            for point_index in 0..npts {
214                let t = unsafe { *self.table.get_unchecked([0, point_index, i, 0]) };
215                for gd in 0..self.gdim {
216                    unsafe {
217                        *points.get_unchecked_mut([gd, point_index]) +=
218                            self.geometry_points.get_value_unchecked([gd, v]) * t
219                    };
220                }
221            }
222        }
223    }
224    fn jacobians<Array3MutImpl: MutableArrayImpl<T, 3>>(
225        &self,
226        entity_index: usize,
227        jacobians: &mut Array<Array3MutImpl, 3>,
228    ) {
229        let npts = self.table.shape()[1];
230        debug_assert!(jacobians.shape()[0] == self.gdim);
231        debug_assert!(jacobians.shape()[1] == self.tdim);
232        debug_assert!(jacobians.shape()[2] == npts);
233
234        jacobians.fill_with_value(T::zero());
235        for i in 0..self.entities.shape()[0] {
236            let v = unsafe { self.entities.get_value_unchecked([i, entity_index]) };
237            for point_index in 0..npts {
238                for td in 0..self.tdim {
239                    let t = unsafe { self.table.get_value_unchecked([1 + td, point_index, i, 0]) };
240                    for gd in 0..self.gdim {
241                        unsafe {
242                            *jacobians.get_unchecked_mut([gd, td, point_index]) +=
243                                self.geometry_points.get_value_unchecked([gd, v]) * t
244                        };
245                    }
246                }
247            }
248        }
249    }
250
251    fn jacobians_inverses_dets<Array3MutImpl: MutableArrayImpl<T, 3>>(
252        &self,
253        entity_index: usize,
254        jacobians: &mut Array<Array3MutImpl, 3>,
255        inverse_jacobians: &mut Array<Array3MutImpl, 3>,
256        jdets: &mut [T],
257    ) {
258        let npts = self.table.shape()[1];
259        debug_assert!(jacobians.shape()[0] == self.gdim);
260        debug_assert!(jacobians.shape()[1] == self.tdim);
261        debug_assert!(jacobians.shape()[2] == npts);
262        debug_assert!(inverse_jacobians.shape()[0] == self.tdim);
263        debug_assert!(inverse_jacobians.shape()[1] == self.gdim);
264        debug_assert!(inverse_jacobians.shape()[2] == npts);
265        debug_assert!(jdets.len() == npts);
266
267        self.jacobians(entity_index, jacobians);
268
269        for point_index in 0..npts {
270            let j = &jacobians.r().slice(2, point_index);
271
272            *jdets.get_mut(point_index).unwrap() =
273                inverse_and_det(j, &mut inverse_jacobians.r_mut().slice(2, point_index));
274        }
275    }
276
277    fn jacobians_inverses_dets_normals<
278        Array2Impl: MutableArrayImpl<T, 2>,
279        Array3MutImpl: MutableArrayImpl<T, 3>,
280    >(
281        &self,
282        entity_index: usize,
283        jacobians: &mut Array<Array3MutImpl, 3>,
284        inverse_jacobians: &mut Array<Array3MutImpl, 3>,
285        jdets: &mut [T],
286        normals: &mut Array<Array2Impl, 2>,
287    ) {
288        if self.tdim + 1 != self.gdim {
289            panic!("Can only compute normal for entities where tdim + 1 == gdim");
290        }
291        let npts = self.table.shape()[1];
292        debug_assert!(jacobians.shape()[0] == self.gdim);
293        debug_assert!(jacobians.shape()[1] == self.tdim);
294        debug_assert!(jacobians.shape()[2] == npts);
295        debug_assert!(inverse_jacobians.shape()[0] == self.tdim);
296        debug_assert!(inverse_jacobians.shape()[1] == self.gdim);
297        debug_assert!(inverse_jacobians.shape()[2] == npts);
298        debug_assert!(jdets.len() == npts);
299        debug_assert!(normals.shape()[0] == self.gdim);
300        debug_assert!(normals.shape()[1] == npts);
301
302        self.jacobians(entity_index, jacobians);
303
304        for point_index in 0..npts {
305            let j = &jacobians.r().slice(2, point_index);
306            let jd = jdets.get_mut(point_index).unwrap();
307
308            *jd = inverse_and_det(j, &mut inverse_jacobians.r_mut().slice(2, point_index));
309
310            let n = &mut normals.r_mut().slice(1, point_index);
311            cross(j, n);
312            for n_i in n.iter_mut() {
313                *n_i /= *jd;
314            }
315        }
316    }
317}
318
319#[cfg(test)]
320mod test {
321    use super::*;
322    use approx::assert_relative_eq;
323    use ndelement::{
324        ciarlet::lagrange,
325        types::{Continuity, ReferenceCellType},
326    };
327    use rand::{Rng, SeedableRng};
328    use rand_chacha::ChaCha8Rng;
329    use rlst::{DynArray, Lu, rlst_dynamic_array};
330
331    fn is_singular(mat: &DynArray<f64, 2>) -> bool {
332        if mat.shape()[0] == 0 || mat.shape()[1] == 0 {
333            return false;
334        }
335        let mut mat_t = rlst_dynamic_array!(f64, [mat.shape()[1], mat.shape()[0]]);
336        mat_t.fill_from(&mat.r().transpose());
337
338        if let Ok(lu) = rlst::dot!(mat_t.r(), mat.r()).lu() {
339            lu.det().abs() < 0.1
340        } else {
341            true
342        }
343    }
344
345    fn non_singular_matrix(gdim: usize, tdim: usize) -> DynArray<f64, 2> {
346        let mut mat = rlst_dynamic_array!(f64, [gdim, tdim]);
347        let mut rng = ChaCha8Rng::seed_from_u64(2);
348        while is_singular(&mat) {
349            for m in mat.iter_mut() {
350                *m = rng.random();
351            }
352        }
353        mat
354    }
355
356    macro_rules! tests {
357        ($gdim:expr, $tdim:expr) => {
358            paste::item! {
359                #[test]
360                fn [< test_inverse_ $gdim _ $tdim >]() {
361                    let mat = non_singular_matrix($gdim, $tdim);
362                    let mut inv = rlst_dynamic_array!(f64, [$tdim, $gdim]);
363
364                    inverse_and_det(&mat, &mut inv);
365
366                    #[allow(clippy::reversed_empty_ranges)]
367                    for i in 0..$tdim {
368                        #[allow(clippy::reversed_empty_ranges)]
369                        for j in 0..$tdim {
370                            let entry = (0..$gdim)
371                                .map(|k| inv[[i, k]] * mat[[k, j]])
372                                .sum::<f64>();
373                            assert_relative_eq!(
374                                entry,
375                                if i == j { 1.0 } else { 0.0 },
376                                epsilon = 1e-10
377                            );
378                        }
379                    }
380
381                }
382
383                #[test]
384                fn [< test_det_ $gdim _ $tdim >]() {
385                    let mat = non_singular_matrix($gdim, $tdim);
386
387                    let rlst_det = if $tdim == 0 || $gdim == 0 {
388                        1.0
389                    } else {
390                        let mut mat_t = rlst_dynamic_array!(f64, [$tdim, $gdim]);
391                        mat_t.fill_from(&mat.r().transpose());
392
393                        if let Ok(lu) = rlst::dot!(mat_t.r(), mat.r()).lu() {
394                            lu.det().sqrt()
395                        } else {
396                            0.0
397                        }
398                    };
399
400                    let mut inv = rlst_dynamic_array!(f64, [$tdim, $gdim]);
401
402                    assert_relative_eq!(
403                        inverse_and_det(&mat, &mut inv),
404                        rlst_det,
405                        epsilon = 1e-10
406                    );
407                }
408            }
409        };
410    }
411
412    tests!(0, 0);
413    tests!(1, 0);
414    tests!(1, 1);
415    tests!(2, 0);
416    tests!(2, 1);
417    tests!(2, 2);
418    tests!(3, 0);
419    tests!(3, 1);
420    tests!(3, 2);
421    tests!(3, 3);
422    tests!(4, 0);
423    tests!(4, 1);
424    tests!(4, 2);
425    tests!(4, 3);
426    tests!(5, 0);
427    tests!(5, 1);
428    tests!(5, 2);
429    tests!(5, 3);
430
431    #[test]
432    fn test_geometry_map_3_2() {
433        let e = lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
434        let mut rng = ChaCha8Rng::seed_from_u64(13);
435
436        let npts = 4;
437        let mut points = rlst_dynamic_array!(f64, [2, npts]);
438
439        for p in 0..npts {
440            *points.get_mut([0, p]).unwrap() = rng.random();
441            *points.get_mut([1, p]).unwrap() =
442                rng.random::<f64>() * points.get_value([0, p]).unwrap();
443        }
444
445        *points.get_mut([0, 0]).unwrap() = 1.0;
446        *points.get_mut([0, 0]).unwrap() = 1.0;
447        let mut geo_points = rlst_dynamic_array!(f64, [3, 3]);
448        *geo_points.get_mut([0, 0]).unwrap() = 1.0;
449        *geo_points.get_mut([1, 0]).unwrap() = 0.0;
450        *geo_points.get_mut([2, 0]).unwrap() = 0.0;
451        *geo_points.get_mut([0, 1]).unwrap() = 2.0;
452        *geo_points.get_mut([1, 1]).unwrap() = 0.0;
453        *geo_points.get_mut([2, 1]).unwrap() = 1.0;
454        *geo_points.get_mut([0, 2]).unwrap() = 0.0;
455        *geo_points.get_mut([1, 2]).unwrap() = 1.0;
456        *geo_points.get_mut([2, 2]).unwrap() = 0.0;
457
458        let mut entities = rlst_dynamic_array!(usize, [2, 1]);
459        *entities.get_mut([0, 0]).unwrap() = 2;
460        *entities.get_mut([1, 0]).unwrap() = 0;
461
462        let gmap = GeometryMap::new(&e, &points, &geo_points, &entities);
463
464        let mut jacobians = rlst_dynamic_array!(f64, [3, 2, npts]);
465        let mut jinv = rlst_dynamic_array!(f64, [2, 3, npts]);
466        let mut jdets = vec![0.0; npts];
467
468        gmap.jacobians_inverses_dets(0, &mut jacobians, &mut jinv, &mut jdets);
469
470        for p in 0..npts {
471            for i in 0..2 {
472                for j in 0..2 {
473                    assert_relative_eq!(
474                        (0..3)
475                            .map(|k| jinv[[i, k, p]] * jacobians[[k, j, p]])
476                            .sum::<f64>(),
477                        if i == j { 1.0 } else { 0.0 },
478                        epsilon = 1e-10
479                    );
480                }
481            }
482        }
483    }
484}
ndgrid/geometry/geometry_map.rs

ndgrid/geometry/
geometry_map.rs
