ndelement/
map.rs

1//! Maps from the reference cell to/from physical cells.
2use crate::traits::Map;
3use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl};
4
5/// Identity map
6///
7/// An identity map pushes values from the reference to the physical
8/// cell without modifying them.
9pub struct IdentityMap {}
10
11impl Map for IdentityMap {
12    fn push_forward<
13        T: RlstScalar,
14        TGeo: RlstScalar,
15        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
16        Array4Impl: ValueArrayImpl<T, 4>,
17        Array4MutImpl: MutableArrayImpl<T, 4>,
18    >(
19        &self,
20        reference_values: &Array<Array4Impl, 4>,
21        nderivs: usize,
22        _jacobians: &Array<Array3GeoImpl, 3>,
23        _jacobian_determinants: &[TGeo],
24        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
25        physical_values: &mut Array<Array4MutImpl, 4>,
26    ) {
27        if nderivs > 0 {
28            unimplemented!();
29        }
30
31        physical_values.fill_from(reference_values);
32    }
33    fn pull_back<
34        T: RlstScalar,
35        TGeo: RlstScalar,
36        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
37        Array4Impl: ValueArrayImpl<T, 4>,
38        Array4MutImpl: MutableArrayImpl<T, 4>,
39    >(
40        &self,
41        physical_values: &Array<Array4Impl, 4>,
42        nderivs: usize,
43        _jacobians: &Array<Array3GeoImpl, 3>,
44        _jacobian_determinants: &[TGeo],
45        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
46        reference_values: &mut Array<Array4MutImpl, 4>,
47    ) {
48        if nderivs > 0 {
49            unimplemented!();
50        }
51
52        reference_values.fill_from(physical_values);
53    }
54    fn physical_value_shape(&self, _gdim: usize) -> Vec<usize> {
55        vec![1]
56    }
57}
58
59/// Covariant Piola map.
60///
61/// Let $F$ be the map from the reference cell to the physical cell
62/// and let $J$ be its Jacobian. Let $\Phi$ be the function values
63/// on the reference cell.  The covariant Piola map is defined by
64/// $$
65/// J^{-T}\Phi\circ F^{-1}
66/// $$
67/// The covariant Piola map preserves tangential continuity. If $J$
68/// is a rectangular matrix then the pseudo-inverse is used instead of
69/// the inverse.
70pub struct CovariantPiolaMap {}
71
72impl Map for CovariantPiolaMap {
73    fn push_forward<
74        T: RlstScalar,
75        TGeo: RlstScalar,
76        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
77        Array4Impl: ValueArrayImpl<T, 4>,
78        Array4MutImpl: MutableArrayImpl<T, 4>,
79    >(
80        &self,
81        reference_values: &Array<Array4Impl, 4>,
82        nderivs: usize,
83        _jacobians: &Array<Array3GeoImpl, 3>,
84        _jacobian_determinants: &[TGeo],
85        inverse_jacobians: &Array<Array3GeoImpl, 3>,
86        physical_values: &mut Array<Array4MutImpl, 4>,
87    ) {
88        let tdim = inverse_jacobians.shape()[0];
89        let gdim = inverse_jacobians.shape()[1];
90        assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
91        assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
92        assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
93        assert_eq!(reference_values.shape()[3], tdim);
94        assert_eq!(physical_values.shape()[3], gdim);
95        if nderivs > 0 {
96            unimplemented!();
97        }
98        for p in 0..reference_values.shape()[1] {
99            for b in 0..reference_values.shape()[2] {
100                for gd in 0..gdim {
101                    unsafe {
102                        *physical_values.get_unchecked_mut([0, p, b, gd]) = (0..tdim)
103                            .map(|td| {
104                                T::from(inverse_jacobians.get_value_unchecked([td, gd, p])).unwrap()
105                                    * reference_values.get_value_unchecked([0, p, b, td])
106                            })
107                            .sum::<T>();
108                    }
109                }
110            }
111        }
112    }
113    fn pull_back<
114        T: RlstScalar,
115        TGeo: RlstScalar,
116        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
117        Array4Impl: ValueArrayImpl<T, 4>,
118        Array4MutImpl: MutableArrayImpl<T, 4>,
119    >(
120        &self,
121        physical_values: &Array<Array4Impl, 4>,
122        nderivs: usize,
123        jacobians: &Array<Array3GeoImpl, 3>,
124        _jacobian_determinants: &[TGeo],
125        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
126        reference_values: &mut Array<Array4MutImpl, 4>,
127    ) {
128        let gdim = jacobians.shape()[0];
129        let tdim = jacobians.shape()[1];
130        assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
131        assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
132        assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
133        assert_eq!(reference_values.shape()[3], tdim);
134        assert_eq!(physical_values.shape()[3], gdim);
135        if nderivs > 0 {
136            unimplemented!();
137        }
138        for p in 0..physical_values.shape()[1] {
139            for b in 0..physical_values.shape()[2] {
140                for td in 0..tdim {
141                    unsafe {
142                        *reference_values.get_unchecked_mut([0, p, b, td]) = (0..gdim)
143                            .map(|gd| {
144                                T::from(jacobians.get_value_unchecked([gd, td, p])).unwrap()
145                                    * physical_values.get_value_unchecked([0, p, b, gd])
146                            })
147                            .sum::<T>();
148                    }
149                }
150            }
151        }
152    }
153    fn physical_value_shape(&self, gdim: usize) -> Vec<usize> {
154        vec![gdim]
155    }
156}
157
158/// Contravariant Piola map.
159///
160/// Let $F$ be the map from the reference cell to the physical cell
161/// and let $J$ be its Jacobian. Let $\Phi$ be the function values
162/// on the reference cell.  The contravariant Piola map is defined by
163/// $$
164/// \frac{1}{\det{J}}J\Phi\circ F^{-1}
165/// $$
166/// The contravariant Piola map preserves normal continuity. If $J$
167/// is a rectangular matrix then $\sqrt{\det{J^TJ}}$ is used instead
168/// of $\det{J}$.
169pub struct ContravariantPiolaMap {}
170
171impl Map for ContravariantPiolaMap {
172    fn push_forward<
173        T: RlstScalar,
174        TGeo: RlstScalar,
175        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
176        Array4Impl: ValueArrayImpl<T, 4>,
177        Array4MutImpl: MutableArrayImpl<T, 4>,
178    >(
179        &self,
180        reference_values: &Array<Array4Impl, 4>,
181        nderivs: usize,
182        jacobians: &Array<Array3GeoImpl, 3>,
183        jacobian_determinants: &[TGeo],
184        _inverse_jacobians: &Array<Array3GeoImpl, 3>,
185        physical_values: &mut Array<Array4MutImpl, 4>,
186    ) {
187        let gdim = jacobians.shape()[0];
188        let tdim = jacobians.shape()[1];
189
190        assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
191        assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
192        assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
193        assert_eq!(reference_values.shape()[3], tdim);
194        assert_eq!(physical_values.shape()[3], gdim);
195        if nderivs > 0 {
196            unimplemented!();
197        }
198        for p in 0..physical_values.shape()[1] {
199            for b in 0..physical_values.shape()[2] {
200                for gd in 0..gdim {
201                    unsafe {
202                        *physical_values.get_unchecked_mut([0, p, b, gd]) = (0..tdim)
203                            .map(|td| {
204                                T::from(jacobians.get_value_unchecked([gd, td, p])).unwrap()
205                                    * reference_values.get_value_unchecked([0, p, b, td])
206                            })
207                            .sum::<T>()
208                            / T::from(*jacobian_determinants.get_unchecked(p)).unwrap();
209                    }
210                }
211            }
212        }
213    }
214    fn pull_back<
215        T: RlstScalar,
216        TGeo: RlstScalar,
217        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
218        Array4Impl: ValueArrayImpl<T, 4>,
219        Array4MutImpl: MutableArrayImpl<T, 4>,
220    >(
221        &self,
222        physical_values: &Array<Array4Impl, 4>,
223        nderivs: usize,
224        _jacobians: &Array<Array3GeoImpl, 3>,
225        jacobian_determinants: &[TGeo],
226        inverse_jacobians: &Array<Array3GeoImpl, 3>,
227        reference_values: &mut Array<Array4MutImpl, 4>,
228    ) {
229        let tdim = inverse_jacobians.shape()[0];
230        let gdim = inverse_jacobians.shape()[1];
231        assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
232        assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
233        assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
234        assert_eq!(reference_values.shape()[3], tdim);
235        assert_eq!(physical_values.shape()[3], gdim);
236        if nderivs > 0 {
237            unimplemented!();
238        }
239        for p in 0..physical_values.shape()[1] {
240            for b in 0..physical_values.shape()[2] {
241                for td in 0..tdim {
242                    unsafe {
243                        *reference_values.get_unchecked_mut([0, p, b, td]) = (0..gdim)
244                            .map(|gd| {
245                                T::from(inverse_jacobians.get_value_unchecked([td, gd, p])).unwrap()
246                                    * physical_values.get_value_unchecked([0, p, b, gd])
247                            })
248                            .sum::<T>()
249                            * T::from(*jacobian_determinants.get_unchecked(p)).unwrap();
250                    }
251                }
252            }
253        }
254    }
255    fn physical_value_shape(&self, gdim: usize) -> Vec<usize> {
256        vec![gdim]
257    }
258}
259
260#[cfg(test)]
261mod test {
262    use super::*;
263    use approx::*;
264    use rlst::{Array, rlst_dynamic_array};
265
266    fn fill_jacobians<T: RlstScalar, Array3MutImpl: MutableArrayImpl<T, 3>>(
267        j: &mut Array<Array3MutImpl, 3>,
268        jdet: &mut [T],
269        jinv: &mut Array<Array3MutImpl, 3>,
270    ) {
271        *j.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
272        *j.get_mut([0, 1, 0]).unwrap() = T::from(1.0).unwrap();
273        *j.get_mut([1, 0, 0]).unwrap() = T::from(0.0).unwrap();
274        *j.get_mut([1, 1, 0]).unwrap() = T::from(1.0).unwrap();
275        *j.get_mut([0, 0, 1]).unwrap() = T::from(2.0).unwrap();
276        *j.get_mut([0, 1, 1]).unwrap() = T::from(0.0).unwrap();
277        *j.get_mut([1, 0, 1]).unwrap() = T::from(0.0).unwrap();
278        *j.get_mut([1, 1, 1]).unwrap() = T::from(3.0).unwrap();
279        jdet[0] = T::from(1.0).unwrap();
280        jdet[1] = T::from(6.0).unwrap();
281        *jinv.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
282        *jinv.get_mut([0, 1, 0]).unwrap() = T::from(-1.0).unwrap();
283        *jinv.get_mut([1, 0, 0]).unwrap() = T::from(0.0).unwrap();
284        *jinv.get_mut([1, 1, 0]).unwrap() = T::from(1.0).unwrap();
285        *jinv.get_mut([0, 0, 1]).unwrap() = T::from(0.5).unwrap();
286        *jinv.get_mut([0, 1, 1]).unwrap() = T::from(0.0).unwrap();
287        *jinv.get_mut([1, 0, 1]).unwrap() = T::from(0.0).unwrap();
288        *jinv.get_mut([1, 1, 1]).unwrap() = T::from(1.0 / 3.0).unwrap();
289    }
290
291    #[test]
292    fn test_identity() {
293        let map = IdentityMap {};
294        let mut values = rlst_dynamic_array!(f64, [1, 2, 2, 1]);
295        *values.get_mut([0, 0, 0, 0]).unwrap() = 1.0;
296        *values.get_mut([0, 1, 0, 0]).unwrap() = 0.0;
297        *values.get_mut([0, 0, 1, 0]).unwrap() = 0.5;
298        *values.get_mut([0, 1, 1, 0]).unwrap() = 2.0;
299
300        let mut j = rlst_dynamic_array!(f64, [2, 2, 2]);
301        let mut jdet = vec![0.0; 2];
302        let mut jinv = rlst_dynamic_array!(f64, [2, 2, 2]);
303        fill_jacobians(&mut j, &mut jdet, &mut jinv);
304
305        let mut physical_values = rlst_dynamic_array!(f64, [1, 2, 2, 1]);
306
307        map.push_forward(&values, 0, &j, &jdet, &jinv, &mut physical_values);
308
309        assert_relative_eq!(
310            *physical_values.get([0, 0, 0, 0]).unwrap(),
311            1.0,
312            epsilon = 1e-14
313        );
314        assert_relative_eq!(
315            *physical_values.get([0, 1, 0, 0]).unwrap(),
316            0.0,
317            epsilon = 1e-14
318        );
319        assert_relative_eq!(
320            *physical_values.get([0, 0, 1, 0]).unwrap(),
321            0.5,
322            epsilon = 1e-14
323        );
324        assert_relative_eq!(
325            *physical_values.get([0, 1, 1, 0]).unwrap(),
326            2.0,
327            epsilon = 1e-14
328        );
329
330        values.set_zero();
331
332        map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values);
333
334        assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14);
335        assert_relative_eq!(*values.get([0, 1, 0, 0]).unwrap(), 0.0, epsilon = 1e-14);
336        assert_relative_eq!(*values.get([0, 0, 1, 0]).unwrap(), 0.5, epsilon = 1e-14);
337        assert_relative_eq!(*values.get([0, 1, 1, 0]).unwrap(), 2.0, epsilon = 1e-14);
338    }
339
340    #[test]
341    fn test_covariant_piola() {
342        let map = CovariantPiolaMap {};
343        let mut values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
344        *values.get_mut([0, 0, 0, 0]).unwrap() = 1.0;
345        *values.get_mut([0, 0, 0, 1]).unwrap() = 0.0;
346        *values.get_mut([0, 1, 0, 0]).unwrap() = 0.0;
347        *values.get_mut([0, 1, 0, 1]).unwrap() = 1.0;
348        *values.get_mut([0, 0, 1, 0]).unwrap() = 0.5;
349        *values.get_mut([0, 0, 1, 1]).unwrap() = 1.5;
350        *values.get_mut([0, 1, 1, 0]).unwrap() = 2.0;
351        *values.get_mut([0, 1, 1, 1]).unwrap() = 2.0;
352
353        let mut j = rlst_dynamic_array!(f64, [2, 2, 2]);
354        let mut jdet = vec![0.0; 2];
355        let mut jinv = rlst_dynamic_array!(f64, [2, 2, 2]);
356        fill_jacobians(&mut j, &mut jdet, &mut jinv);
357
358        let mut physical_values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
359
360        map.push_forward(&values, 0, &j, &jdet, &jinv, &mut physical_values);
361
362        assert_relative_eq!(
363            *physical_values.get([0, 0, 0, 0]).unwrap(),
364            1.0,
365            epsilon = 1e-14
366        );
367        assert_relative_eq!(
368            *physical_values.get([0, 0, 0, 1]).unwrap(),
369            -1.0,
370            epsilon = 1e-14
371        );
372        assert_relative_eq!(
373            *physical_values.get([0, 1, 0, 0]).unwrap(),
374            0.0,
375            epsilon = 1e-14
376        );
377        assert_relative_eq!(
378            *physical_values.get([0, 1, 0, 1]).unwrap(),
379            1.0 / 3.0,
380            epsilon = 1e-14
381        );
382        assert_relative_eq!(
383            *physical_values.get([0, 0, 1, 0]).unwrap(),
384            0.5,
385            epsilon = 1e-14
386        );
387        assert_relative_eq!(
388            *physical_values.get([0, 0, 1, 1]).unwrap(),
389            1.0,
390            epsilon = 1e-14
391        );
392        assert_relative_eq!(
393            *physical_values.get([0, 1, 1, 0]).unwrap(),
394            1.0,
395            epsilon = 1e-14
396        );
397        assert_relative_eq!(
398            *physical_values.get([0, 1, 1, 1]).unwrap(),
399            2.0 / 3.0,
400            epsilon = 1e-14
401        );
402
403        values.set_zero();
404        map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values);
405
406        assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14);
407        assert_relative_eq!(*values.get([0, 0, 0, 1]).unwrap(), 0.0, epsilon = 1e-14);
408        assert_relative_eq!(*values.get([0, 1, 0, 0]).unwrap(), 0.0, epsilon = 1e-14);
409        assert_relative_eq!(*values.get([0, 1, 0, 1]).unwrap(), 1.0, epsilon = 1e-14);
410        assert_relative_eq!(*values.get([0, 0, 1, 0]).unwrap(), 0.5, epsilon = 1e-14);
411        assert_relative_eq!(*values.get([0, 0, 1, 1]).unwrap(), 1.5, epsilon = 1e-14);
412        assert_relative_eq!(*values.get([0, 1, 1, 0]).unwrap(), 2.0, epsilon = 1e-14);
413        assert_relative_eq!(*values.get([0, 1, 1, 1]).unwrap(), 2.0, epsilon = 1e-14);
414    }
415
416    #[test]
417    fn test_contravariant_piola() {
418        let map = ContravariantPiolaMap {};
419        let mut values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
420        *values.get_mut([0, 0, 0, 0]).unwrap() = 1.0;
421        *values.get_mut([0, 0, 0, 1]).unwrap() = 0.0;
422        *values.get_mut([0, 1, 0, 0]).unwrap() = 0.0;
423        *values.get_mut([0, 1, 0, 1]).unwrap() = 1.0;
424        *values.get_mut([0, 0, 1, 0]).unwrap() = 0.5;
425        *values.get_mut([0, 0, 1, 1]).unwrap() = 1.5;
426        *values.get_mut([0, 1, 1, 0]).unwrap() = 2.0;
427        *values.get_mut([0, 1, 1, 1]).unwrap() = 2.0;
428
429        let mut j = rlst_dynamic_array!(f64, [2, 2, 2]);
430        let mut jdet = vec![0.0; 2];
431        let mut jinv = rlst_dynamic_array!(f64, [2, 2, 2]);
432        fill_jacobians(&mut j, &mut jdet, &mut jinv);
433
434        let mut physical_values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
435
436        map.push_forward(&values, 0, &j, &jdet, &jinv, &mut physical_values);
437
438        assert_relative_eq!(
439            *physical_values.get([0, 0, 0, 0]).unwrap(),
440            1.0,
441            epsilon = 1e-14
442        );
443        assert_relative_eq!(
444            *physical_values.get([0, 0, 0, 1]).unwrap(),
445            0.0,
446            epsilon = 1e-14
447        );
448        assert_relative_eq!(
449            *physical_values.get([0, 1, 0, 0]).unwrap(),
450            0.0,
451            epsilon = 1e-14
452        );
453        assert_relative_eq!(
454            *physical_values.get([0, 1, 0, 1]).unwrap(),
455            0.5,
456            epsilon = 1e-14
457        );
458        assert_relative_eq!(
459            *physical_values.get([0, 0, 1, 0]).unwrap(),
460            2.0,
461            epsilon = 1e-14
462        );
463        assert_relative_eq!(
464            *physical_values.get([0, 0, 1, 1]).unwrap(),
465            1.5,
466            epsilon = 1e-14
467        );
468        assert_relative_eq!(
469            *physical_values.get([0, 1, 1, 0]).unwrap(),
470            2.0 / 3.0,
471            epsilon = 1e-14
472        );
473        assert_relative_eq!(
474            *physical_values.get([0, 1, 1, 1]).unwrap(),
475            1.0,
476            epsilon = 1e-14
477        );
478
479        values.set_zero();
480        map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values);
481
482        assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14);
483        assert_relative_eq!(*values.get([0, 0, 0, 1]).unwrap(), 0.0, epsilon = 1e-14);
484        assert_relative_eq!(*values.get([0, 1, 0, 0]).unwrap(), 0.0, epsilon = 1e-14);
485        assert_relative_eq!(*values.get([0, 1, 0, 1]).unwrap(), 1.0, epsilon = 1e-14);
486        assert_relative_eq!(*values.get([0, 0, 1, 0]).unwrap(), 0.5, epsilon = 1e-14);
487        assert_relative_eq!(*values.get([0, 0, 1, 1]).unwrap(), 1.5, epsilon = 1e-14);
488        assert_relative_eq!(*values.get([0, 1, 1, 0]).unwrap(), 2.0, epsilon = 1e-14);
489        assert_relative_eq!(*values.get([0, 1, 1, 1]).unwrap(), 2.0, epsilon = 1e-14);
490    }
491}
ndelement/map.rs

ndelement/
map.rs
