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()[1];
89        let gdim = inverse_jacobians.shape()[2];
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([p, td, gd])).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()[1];
129        let tdim = jacobians.shape()[2];
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([p, gd, td])).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()[1];
188        let tdim = jacobians.shape()[2];
189        assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
190        assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
191        assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
192        assert_eq!(reference_values.shape()[3], tdim);
193        assert_eq!(physical_values.shape()[3], gdim);
194        if nderivs > 0 {
195            unimplemented!();
196        }
197        for p in 0..physical_values.shape()[1] {
198            for b in 0..physical_values.shape()[2] {
199                for gd in 0..gdim {
200                    unsafe {
201                        *physical_values.get_unchecked_mut([0, p, b, gd]) = (0..tdim)
202                            .map(|td| {
203                                T::from(jacobians.get_value_unchecked([p, gd, td])).unwrap()
204                                    * reference_values.get_value_unchecked([0, p, b, td])
205                            })
206                            .sum::<T>()
207                            / T::from(*jacobian_determinants.get_unchecked(p)).unwrap();
208                    }
209                }
210            }
211        }
212    }
213    fn pull_back<
214        T: RlstScalar,
215        TGeo: RlstScalar,
216        Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
217        Array4Impl: ValueArrayImpl<T, 4>,
218        Array4MutImpl: MutableArrayImpl<T, 4>,
219    >(
220        &self,
221        physical_values: &Array<Array4Impl, 4>,
222        nderivs: usize,
223        _jacobians: &Array<Array3GeoImpl, 3>,
224        jacobian_determinants: &[TGeo],
225        inverse_jacobians: &Array<Array3GeoImpl, 3>,
226        reference_values: &mut Array<Array4MutImpl, 4>,
227    ) {
228        let tdim = inverse_jacobians.shape()[1];
229        let gdim = inverse_jacobians.shape()[2];
230        assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
231        assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
232        assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
233        assert_eq!(reference_values.shape()[3], tdim);
234        assert_eq!(physical_values.shape()[3], gdim);
235        if nderivs > 0 {
236            unimplemented!();
237        }
238        for p in 0..physical_values.shape()[1] {
239            for b in 0..physical_values.shape()[2] {
240                for td in 0..tdim {
241                    unsafe {
242                        *reference_values.get_unchecked_mut([0, p, b, td]) = (0..gdim)
243                            .map(|gd| {
244                                T::from(inverse_jacobians.get_value_unchecked([p, td, gd])).unwrap()
245                                    * physical_values.get_value_unchecked([0, p, b, gd])
246                            })
247                            .sum::<T>()
248                            * T::from(*jacobian_determinants.get_unchecked(p)).unwrap();
249                    }
250                }
251            }
252        }
253    }
254    fn physical_value_shape(&self, gdim: usize) -> Vec<usize> {
255        vec![gdim]
256    }
257}
258
259#[cfg(test)]
260mod test {
261    use super::*;
262    use approx::*;
263    use rlst::{Array, rlst_dynamic_array};
264
265    fn fill_jacobians<T: RlstScalar, Array3MutImpl: MutableArrayImpl<T, 3>>(
266        j: &mut Array<Array3MutImpl, 3>,
267        jdet: &mut [T],
268        jinv: &mut Array<Array3MutImpl, 3>,
269    ) {
270        *j.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
271        *j.get_mut([0, 0, 1]).unwrap() = T::from(1.0).unwrap();
272        *j.get_mut([0, 1, 0]).unwrap() = T::from(0.0).unwrap();
273        *j.get_mut([0, 1, 1]).unwrap() = T::from(1.0).unwrap();
274        *j.get_mut([1, 0, 0]).unwrap() = T::from(2.0).unwrap();
275        *j.get_mut([1, 0, 1]).unwrap() = T::from(0.0).unwrap();
276        *j.get_mut([1, 1, 0]).unwrap() = T::from(0.0).unwrap();
277        *j.get_mut([1, 1, 1]).unwrap() = T::from(3.0).unwrap();
278        jdet[0] = T::from(1.0).unwrap();
279        jdet[1] = T::from(6.0).unwrap();
280        *jinv.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
281        *jinv.get_mut([0, 0, 1]).unwrap() = T::from(-1.0).unwrap();
282        *jinv.get_mut([0, 1, 0]).unwrap() = T::from(0.0).unwrap();
283        *jinv.get_mut([0, 1, 1]).unwrap() = T::from(1.0).unwrap();
284        *jinv.get_mut([1, 0, 0]).unwrap() = T::from(0.5).unwrap();
285        *jinv.get_mut([1, 0, 1]).unwrap() = T::from(0.0).unwrap();
286        *jinv.get_mut([1, 1, 0]).unwrap() = T::from(0.0).unwrap();
287        *jinv.get_mut([1, 1, 1]).unwrap() = T::from(1.0 / 3.0).unwrap();
288    }
289
290    #[test]
291    fn test_identity() {
292        let map = IdentityMap {};
293        let mut values = rlst_dynamic_array!(f64, [1, 2, 2, 1]);
294        *values.get_mut([0, 0, 0, 0]).unwrap() = 1.0;
295        *values.get_mut([0, 1, 0, 0]).unwrap() = 0.0;
296        *values.get_mut([0, 0, 1, 0]).unwrap() = 0.5;
297        *values.get_mut([0, 1, 1, 0]).unwrap() = 2.0;
298
299        let mut j = rlst_dynamic_array!(f64, [2, 2, 2]);
300        let mut jdet = vec![0.0; 2];
301        let mut jinv = rlst_dynamic_array!(f64, [2, 2, 2]);
302        fill_jacobians(&mut j, &mut jdet, &mut jinv);
303
304        let mut physical_values = rlst_dynamic_array!(f64, [1, 2, 2, 1]);
305
306        map.push_forward(&values, 0, &j, &jdet, &jinv, &mut physical_values);
307
308        assert_relative_eq!(
309            *physical_values.get([0, 0, 0, 0]).unwrap(),
310            1.0,
311            epsilon = 1e-14
312        );
313        assert_relative_eq!(
314            *physical_values.get([0, 1, 0, 0]).unwrap(),
315            0.0,
316            epsilon = 1e-14
317        );
318        assert_relative_eq!(
319            *physical_values.get([0, 0, 1, 0]).unwrap(),
320            0.5,
321            epsilon = 1e-14
322        );
323        assert_relative_eq!(
324            *physical_values.get([0, 1, 1, 0]).unwrap(),
325            2.0,
326            epsilon = 1e-14
327        );
328
329        values.set_zero();
330
331        map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values);
332
333        assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14);
334        assert_relative_eq!(*values.get([0, 1, 0, 0]).unwrap(), 0.0, epsilon = 1e-14);
335        assert_relative_eq!(*values.get([0, 0, 1, 0]).unwrap(), 0.5, epsilon = 1e-14);
336        assert_relative_eq!(*values.get([0, 1, 1, 0]).unwrap(), 2.0, epsilon = 1e-14);
337    }
338
339    #[test]
340    fn test_covariant_piola() {
341        let map = CovariantPiolaMap {};
342        let mut values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
343        *values.get_mut([0, 0, 0, 0]).unwrap() = 1.0;
344        *values.get_mut([0, 0, 0, 1]).unwrap() = 0.0;
345        *values.get_mut([0, 1, 0, 0]).unwrap() = 0.0;
346        *values.get_mut([0, 1, 0, 1]).unwrap() = 1.0;
347        *values.get_mut([0, 0, 1, 0]).unwrap() = 0.5;
348        *values.get_mut([0, 0, 1, 1]).unwrap() = 1.5;
349        *values.get_mut([0, 1, 1, 0]).unwrap() = 2.0;
350        *values.get_mut([0, 1, 1, 1]).unwrap() = 2.0;
351
352        let mut j = rlst_dynamic_array!(f64, [2, 2, 2]);
353        let mut jdet = vec![0.0; 2];
354        let mut jinv = rlst_dynamic_array!(f64, [2, 2, 2]);
355        fill_jacobians(&mut j, &mut jdet, &mut jinv);
356
357        let mut physical_values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
358
359        map.push_forward(&values, 0, &j, &jdet, &jinv, &mut physical_values);
360
361        assert_relative_eq!(
362            *physical_values.get([0, 0, 0, 0]).unwrap(),
363            1.0,
364            epsilon = 1e-14
365        );
366        assert_relative_eq!(
367            *physical_values.get([0, 0, 0, 1]).unwrap(),
368            -1.0,
369            epsilon = 1e-14
370        );
371        assert_relative_eq!(
372            *physical_values.get([0, 1, 0, 0]).unwrap(),
373            0.0,
374            epsilon = 1e-14
375        );
376        assert_relative_eq!(
377            *physical_values.get([0, 1, 0, 1]).unwrap(),
378            1.0 / 3.0,
379            epsilon = 1e-14
380        );
381        assert_relative_eq!(
382            *physical_values.get([0, 0, 1, 0]).unwrap(),
383            0.5,
384            epsilon = 1e-14
385        );
386        assert_relative_eq!(
387            *physical_values.get([0, 0, 1, 1]).unwrap(),
388            1.0,
389            epsilon = 1e-14
390        );
391        assert_relative_eq!(
392            *physical_values.get([0, 1, 1, 0]).unwrap(),
393            1.0,
394            epsilon = 1e-14
395        );
396        assert_relative_eq!(
397            *physical_values.get([0, 1, 1, 1]).unwrap(),
398            2.0 / 3.0,
399            epsilon = 1e-14
400        );
401
402        values.set_zero();
403        map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values);
404
405        assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14);
406        assert_relative_eq!(*values.get([0, 0, 0, 1]).unwrap(), 0.0, epsilon = 1e-14);
407        assert_relative_eq!(*values.get([0, 1, 0, 0]).unwrap(), 0.0, epsilon = 1e-14);
408        assert_relative_eq!(*values.get([0, 1, 0, 1]).unwrap(), 1.0, epsilon = 1e-14);
409        assert_relative_eq!(*values.get([0, 0, 1, 0]).unwrap(), 0.5, epsilon = 1e-14);
410        assert_relative_eq!(*values.get([0, 0, 1, 1]).unwrap(), 1.5, epsilon = 1e-14);
411        assert_relative_eq!(*values.get([0, 1, 1, 0]).unwrap(), 2.0, epsilon = 1e-14);
412        assert_relative_eq!(*values.get([0, 1, 1, 1]).unwrap(), 2.0, epsilon = 1e-14);
413    }
414
415    #[test]
416    fn test_contravariant_piola() {
417        let map = ContravariantPiolaMap {};
418        let mut values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
419        *values.get_mut([0, 0, 0, 0]).unwrap() = 1.0;
420        *values.get_mut([0, 0, 0, 1]).unwrap() = 0.0;
421        *values.get_mut([0, 1, 0, 0]).unwrap() = 0.0;
422        *values.get_mut([0, 1, 0, 1]).unwrap() = 1.0;
423        *values.get_mut([0, 0, 1, 0]).unwrap() = 0.5;
424        *values.get_mut([0, 0, 1, 1]).unwrap() = 1.5;
425        *values.get_mut([0, 1, 1, 0]).unwrap() = 2.0;
426        *values.get_mut([0, 1, 1, 1]).unwrap() = 2.0;
427
428        let mut j = rlst_dynamic_array!(f64, [2, 2, 2]);
429        let mut jdet = vec![0.0; 2];
430        let mut jinv = rlst_dynamic_array!(f64, [2, 2, 2]);
431        fill_jacobians(&mut j, &mut jdet, &mut jinv);
432
433        let mut physical_values = rlst_dynamic_array!(f64, [1, 2, 2, 2]);
434
435        map.push_forward(&values, 0, &j, &jdet, &jinv, &mut physical_values);
436
437        assert_relative_eq!(
438            *physical_values.get([0, 0, 0, 0]).unwrap(),
439            1.0,
440            epsilon = 1e-14
441        );
442        assert_relative_eq!(
443            *physical_values.get([0, 0, 0, 1]).unwrap(),
444            0.0,
445            epsilon = 1e-14
446        );
447        assert_relative_eq!(
448            *physical_values.get([0, 1, 0, 0]).unwrap(),
449            0.0,
450            epsilon = 1e-14
451        );
452        assert_relative_eq!(
453            *physical_values.get([0, 1, 0, 1]).unwrap(),
454            0.5,
455            epsilon = 1e-14
456        );
457        assert_relative_eq!(
458            *physical_values.get([0, 0, 1, 0]).unwrap(),
459            2.0,
460            epsilon = 1e-14
461        );
462        assert_relative_eq!(
463            *physical_values.get([0, 0, 1, 1]).unwrap(),
464            1.5,
465            epsilon = 1e-14
466        );
467        assert_relative_eq!(
468            *physical_values.get([0, 1, 1, 0]).unwrap(),
469            2.0 / 3.0,
470            epsilon = 1e-14
471        );
472        assert_relative_eq!(
473            *physical_values.get([0, 1, 1, 1]).unwrap(),
474            1.0,
475            epsilon = 1e-14
476        );
477
478        values.set_zero();
479        map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values);
480
481        assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14);
482        assert_relative_eq!(*values.get([0, 0, 0, 1]).unwrap(), 0.0, epsilon = 1e-14);
483        assert_relative_eq!(*values.get([0, 1, 0, 0]).unwrap(), 0.0, epsilon = 1e-14);
484        assert_relative_eq!(*values.get([0, 1, 0, 1]).unwrap(), 1.0, epsilon = 1e-14);
485        assert_relative_eq!(*values.get([0, 0, 1, 0]).unwrap(), 0.5, epsilon = 1e-14);
486        assert_relative_eq!(*values.get([0, 0, 1, 1]).unwrap(), 1.5, epsilon = 1e-14);
487        assert_relative_eq!(*values.get([0, 1, 1, 0]).unwrap(), 2.0, epsilon = 1e-14);
488        assert_relative_eq!(*values.get([0, 1, 1, 1]).unwrap(), 2.0, epsilon = 1e-14);
489    }
490}
ndelement/map.rs

ndelement/
map.rs
