ndelement/ciarlet/
lagrange.rs

1//! Lagrange elements.
2
3use super::CiarletElement;
4use crate::map::IdentityMap;
5use crate::polynomials::polynomial_count;
6use crate::reference_cell;
7use crate::traits::ElementFamily;
8use crate::types::{Continuity, ReferenceCellType};
9use quadraturerules::{Domain, QuadratureRule, single_integral_quadrature};
10use rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri};
11use rlst::{RlstScalar, rlst_dynamic_array};
12use std::marker::PhantomData;
13
14/// Variant of a Lagrange element
15#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
16#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)]
17#[repr(u8)]
18pub enum LagrangeVariant {
19    /// Element defined using equispaced points
20    Equispaced,
21    /// Element defined using Guass-Lobatto-Legendre (GLL) points
22    GLL,
23}
24
25/// Create a Lagrange element.
26pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
27    cell_type: ReferenceCellType,
28    degree: usize,
29    continuity: Continuity,
30    variant: LagrangeVariant,
31) -> CiarletElement<T, IdentityMap, TGeo> {
32    let dim = polynomial_count(cell_type, degree);
33    let tdim = reference_cell::dim(cell_type);
34    let mut wcoeffs = rlst_dynamic_array!(T, [dim, 1, dim]);
35    for i in 0..dim {
36        *wcoeffs.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap();
37    }
38
39    let mut x = [vec![], vec![], vec![], vec![]];
40    let mut m = [vec![], vec![], vec![], vec![]];
41    let entity_counts = reference_cell::entity_counts(cell_type);
42    let vertices = reference_cell::vertices::<TGeo>(cell_type);
43    if degree == 0 {
44        if continuity == Continuity::Standard {
45            panic!("Cannot create continuous degree 0 Lagrange element");
46        }
47        for (d, counts) in entity_counts.iter().enumerate() {
48            for _e in 0..*counts {
49                x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
50                m[d].push(rlst_dynamic_array!(T, [0, 1, 0]));
51            }
52        }
53        let mut midp = rlst_dynamic_array!(TGeo, [tdim, 1]);
54        let nvertices = entity_counts[0];
55        for i in 0..tdim {
56            for vertex in &vertices {
57                *midp.get_mut([i, 0]).unwrap() += num::cast::<_, TGeo>(vertex[i]).unwrap();
58            }
59            *midp.get_mut([i, 0]).unwrap() /= num::cast::<_, TGeo>(nvertices).unwrap();
60        }
61        x[tdim].push(midp);
62        let mut mentry = rlst_dynamic_array!(T, [1, 1, 1]);
63        *mentry.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
64        m[tdim].push(mentry);
65    } else {
66        let edges = reference_cell::edges(cell_type);
67        let faces = reference_cell::faces(cell_type);
68        let volumes = reference_cell::volumes(cell_type);
69        for vertex in &vertices {
70            let mut pts = rlst_dynamic_array!(TGeo, [tdim, 1]);
71            for (i, v) in vertex.iter().enumerate() {
72                *pts.get_mut([i, 0]).unwrap() = num::cast::<_, TGeo>(*v).unwrap();
73            }
74            x[0].push(pts);
75            let mut mentry = rlst_dynamic_array!(T, [1, 1, 1]);
76            *mentry.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
77            m[0].push(mentry);
78        }
79
80        let pts1d = match variant {
81            LagrangeVariant::Equispaced => (1..degree)
82                .map(|i| num::cast::<_, TGeo>(i).unwrap() / num::cast::<_, TGeo>(degree).unwrap())
83                .collect::<Vec<_>>(),
84            LagrangeVariant::GLL => {
85                let (points, _weights) = single_integral_quadrature(
86                    QuadratureRule::GaussLobattoLegendre,
87                    Domain::Interval,
88                    degree - 1,
89                )
90                .unwrap();
91                (1..degree)
92                    .map(|i| num::cast::<_, TGeo>(points[2 * i + 1]).unwrap())
93                    .collect::<Vec<_>>()
94            }
95        };
96
97        for e in &edges {
98            let mut pts = rlst_dynamic_array!(TGeo, [tdim, degree - 1]);
99            let [vn0, vn1] = e[..] else {
100                panic!();
101            };
102            let v0 = &vertices[vn0];
103            let v1 = &vertices[vn1];
104            let mut ident = rlst_dynamic_array!(T, [degree - 1, 1, degree - 1]);
105
106            for (i, p) in pts1d.iter().enumerate() {
107                *ident.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap();
108                for j in 0..tdim {
109                    *pts.get_mut([j, i]).unwrap() = num::cast::<_, TGeo>(v0[j]).unwrap()
110                        + *p * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap();
111                }
112            }
113            x[1].push(pts);
114            m[1].push(ident);
115        }
116        for (e, face_type) in reference_cell::entity_types(cell_type)[2]
117            .iter()
118            .enumerate()
119        {
120            let npts = match face_type {
121                ReferenceCellType::Triangle => {
122                    if degree > 2 {
123                        (degree - 1) * (degree - 2) / 2
124                    } else {
125                        0
126                    }
127                }
128                ReferenceCellType::Quadrilateral => (degree - 1).pow(2),
129                _ => {
130                    panic!("Unsupported face type");
131                }
132            };
133            let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]);
134
135            let [vn0, vn1, vn2] = faces[e][..3] else {
136                panic!();
137            };
138            let v0 = &vertices[vn0];
139            let v1 = &vertices[vn1];
140            let v2 = &vertices[vn2];
141
142            match face_type {
143                ReferenceCellType::Triangle => {
144                    if variant != LagrangeVariant::Equispaced {
145                        unimplemented!();
146                    }
147                    let mut n = 0;
148                    for i1 in 1..degree {
149                        for i0 in 1..degree - i1 {
150                            for j in 0..tdim {
151                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
152                                    .unwrap()
153                                    + num::cast::<_, TGeo>(i0).unwrap()
154                                        / num::cast::<_, TGeo>(degree).unwrap()
155                                        * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
156                                    + num::cast::<_, TGeo>(i1).unwrap()
157                                        / num::cast::<_, TGeo>(degree).unwrap()
158                                        * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap();
159                            }
160                            n += 1;
161                        }
162                    }
163                }
164                ReferenceCellType::Quadrilateral => {
165                    let mut n = 0;
166                    for p1 in &pts1d {
167                        for p0 in &pts1d {
168                            for j in 0..tdim {
169                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
170                                    .unwrap()
171                                    + *p0 * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
172                                    + *p1 * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap();
173                            }
174                            n += 1;
175                        }
176                    }
177                }
178                _ => {
179                    panic!("Unsupported face type.");
180                }
181            };
182
183            let mut ident = rlst_dynamic_array!(T, [npts, 1, npts]);
184            for i in 0..npts {
185                *ident.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap();
186            }
187            x[2].push(pts);
188            m[2].push(ident);
189        }
190        for (e, volume_type) in reference_cell::entity_types(cell_type)[3]
191            .iter()
192            .enumerate()
193        {
194            let npts = match volume_type {
195                ReferenceCellType::Tetrahedron => {
196                    if degree > 2 {
197                        (degree - 1) * (degree - 2) * (degree - 3) / 6
198                    } else {
199                        0
200                    }
201                }
202                ReferenceCellType::Hexahedron => (degree - 1).pow(3),
203                _ => {
204                    panic!("Unsupported face type");
205                }
206            };
207            let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]);
208
209            match volume_type {
210                ReferenceCellType::Tetrahedron => {
211                    if variant != LagrangeVariant::Equispaced {
212                        unimplemented!();
213                    }
214                    let [vn0, vn1, vn2, vn3] = volumes[e][..4] else {
215                        panic!();
216                    };
217                    let v0 = &vertices[vn0];
218                    let v1 = &vertices[vn1];
219                    let v2 = &vertices[vn2];
220                    let v3 = &vertices[vn3];
221
222                    let mut n = 0;
223                    for i2 in 1..degree {
224                        for i1 in 1..degree - i2 {
225                            for i0 in 1..degree - i2 - i1 {
226                                for j in 0..tdim {
227                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
228                                        .unwrap()
229                                        + num::cast::<_, TGeo>(i0).unwrap()
230                                            / num::cast::<_, TGeo>(degree).unwrap()
231                                            * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
232                                        + num::cast::<_, TGeo>(i1).unwrap()
233                                            / num::cast::<_, TGeo>(degree).unwrap()
234                                            * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap()
235                                        + num::cast::<_, TGeo>(i2).unwrap()
236                                            / num::cast::<_, TGeo>(degree).unwrap()
237                                            * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap();
238                                }
239                                n += 1;
240                            }
241                        }
242                    }
243                }
244                ReferenceCellType::Hexahedron => {
245                    let [vn0, vn1, vn2, _, vn3] = volumes[e][..5] else {
246                        panic!();
247                    };
248                    let v0 = &vertices[vn0];
249                    let v1 = &vertices[vn1];
250                    let v2 = &vertices[vn2];
251                    let v3 = &vertices[vn3];
252
253                    let mut n = 0;
254                    for p2 in &pts1d {
255                        for p1 in &pts1d {
256                            for p0 in &pts1d {
257                                for j in 0..tdim {
258                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
259                                        .unwrap()
260                                        + *p0 * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
261                                        + *p1 * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap()
262                                        + *p2 * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap();
263                                }
264                                n += 1;
265                            }
266                        }
267                    }
268                }
269                _ => {
270                    panic!("Unsupported face type.");
271                }
272            };
273
274            let mut ident = rlst_dynamic_array!(T, [npts, 1, npts]);
275            for i in 0..npts {
276                *ident.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap();
277            }
278            x[3].push(pts);
279            m[3].push(ident);
280        }
281    }
282    CiarletElement::<T, IdentityMap, TGeo>::create(
283        "Lagrange".to_string(),
284        cell_type,
285        degree,
286        vec![],
287        wcoeffs,
288        x,
289        m,
290        continuity,
291        degree,
292        IdentityMap {},
293    )
294}
295
296/// Lagrange element family.
297///
298/// A family of Lagrange elements on multiple cell types with appropriate
299/// continuity across different cell types.
300pub struct LagrangeElementFamily<T: RlstScalar + Getrf + Getri = f64, TGeo: RlstScalar = f64> {
301    degree: usize,
302    continuity: Continuity,
303    variant: LagrangeVariant,
304    _t: PhantomData<T>,
305    _tgeo: PhantomData<TGeo>,
306}
307
308impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> LagrangeElementFamily<T, TGeo> {
309    /// Create new family with given `degree` and `continuity`.
310    pub fn new(degree: usize, continuity: Continuity, variant: LagrangeVariant) -> Self {
311        Self {
312            degree,
313            continuity,
314            variant,
315            _t: PhantomData,
316            _tgeo: PhantomData,
317        }
318    }
319}
320
321impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> ElementFamily
322    for LagrangeElementFamily<T, TGeo>
323{
324    type T = T;
325    type FiniteElement = CiarletElement<T, IdentityMap, TGeo>;
326    type CellType = ReferenceCellType;
327    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, IdentityMap, TGeo> {
328        create::<T, TGeo>(cell_type, self.degree, self.continuity, self.variant)
329    }
330}
331
332#[cfg(test)]
333mod test {
334    use super::super::test::check_dofs;
335    use super::*;
336    use crate::traits::FiniteElement;
337    use approx::*;
338    use paste::paste;
339    use rlst::DynArray;
340
341    #[test]
342    fn test_lagrange_1() {
343        let e = create::<f64, f64>(
344            ReferenceCellType::Triangle,
345            1,
346            Continuity::Standard,
347            LagrangeVariant::Equispaced,
348        );
349        assert_eq!(e.value_size(), 1);
350    }
351
352    #[test]
353    fn test_lagrange_0_interval() {
354        let e = create::<f64, f64>(
355            ReferenceCellType::Interval,
356            0,
357            Continuity::Discontinuous,
358            LagrangeVariant::Equispaced,
359        );
360        assert_eq!(e.value_size(), 1);
361        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
362        let mut points = rlst_dynamic_array!(f64, [1, 4]);
363        *points.get_mut([0, 0]).unwrap() = 0.0;
364        *points.get_mut([0, 1]).unwrap() = 0.2;
365        *points.get_mut([0, 2]).unwrap() = 0.4;
366        *points.get_mut([0, 3]).unwrap() = 1.0;
367        e.tabulate(&points, 0, &mut data);
368
369        for pt in 0..4 {
370            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
371        }
372        check_dofs(e);
373    }
374
375    #[test]
376    fn test_lagrange_1_interval() {
377        let e = create::<f64, f64>(
378            ReferenceCellType::Interval,
379            1,
380            Continuity::Standard,
381            LagrangeVariant::Equispaced,
382        );
383        assert_eq!(e.value_size(), 1);
384        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
385        let mut points = rlst_dynamic_array!(f64, [1, 4]);
386        *points.get_mut([0, 0]).unwrap() = 0.0;
387        *points.get_mut([0, 1]).unwrap() = 0.2;
388        *points.get_mut([0, 2]).unwrap() = 0.4;
389        *points.get_mut([0, 3]).unwrap() = 1.0;
390        e.tabulate(&points, 0, &mut data);
391
392        for pt in 0..4 {
393            let x = *points.get([0, pt]).unwrap();
394            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0 - x);
395            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x);
396        }
397        check_dofs(e);
398    }
399
400    #[test]
401    fn test_lagrange_0_triangle() {
402        let e = create::<f64, f64>(
403            ReferenceCellType::Triangle,
404            0,
405            Continuity::Discontinuous,
406            LagrangeVariant::Equispaced,
407        );
408        assert_eq!(e.value_size(), 1);
409        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
410
411        let mut points = rlst_dynamic_array!(f64, [2, 6]);
412        *points.get_mut([0, 0]).unwrap() = 0.0;
413        *points.get_mut([1, 0]).unwrap() = 0.0;
414        *points.get_mut([0, 1]).unwrap() = 1.0;
415        *points.get_mut([1, 1]).unwrap() = 0.0;
416        *points.get_mut([0, 2]).unwrap() = 0.0;
417        *points.get_mut([1, 2]).unwrap() = 1.0;
418        *points.get_mut([0, 3]).unwrap() = 0.5;
419        *points.get_mut([1, 3]).unwrap() = 0.0;
420        *points.get_mut([0, 4]).unwrap() = 0.0;
421        *points.get_mut([1, 4]).unwrap() = 0.5;
422        *points.get_mut([0, 5]).unwrap() = 0.5;
423        *points.get_mut([1, 5]).unwrap() = 0.5;
424
425        e.tabulate(&points, 0, &mut data);
426
427        for pt in 0..6 {
428            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
429        }
430        check_dofs(e);
431    }
432
433    #[test]
434    fn test_lagrange_1_triangle() {
435        let e = create::<f64, f64>(
436            ReferenceCellType::Triangle,
437            1,
438            Continuity::Standard,
439            LagrangeVariant::Equispaced,
440        );
441        assert_eq!(e.value_size(), 1);
442        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
443        let mut points = rlst_dynamic_array!(f64, [2, 6]);
444        *points.get_mut([0, 0]).unwrap() = 0.0;
445        *points.get_mut([1, 0]).unwrap() = 0.0;
446        *points.get_mut([0, 1]).unwrap() = 1.0;
447        *points.get_mut([1, 1]).unwrap() = 0.0;
448        *points.get_mut([0, 2]).unwrap() = 0.0;
449        *points.get_mut([1, 2]).unwrap() = 1.0;
450        *points.get_mut([0, 3]).unwrap() = 0.5;
451        *points.get_mut([1, 3]).unwrap() = 0.0;
452        *points.get_mut([0, 4]).unwrap() = 0.0;
453        *points.get_mut([1, 4]).unwrap() = 0.5;
454        *points.get_mut([0, 5]).unwrap() = 0.5;
455        *points.get_mut([1, 5]).unwrap() = 0.5;
456        e.tabulate(&points, 0, &mut data);
457
458        for pt in 0..6 {
459            let x = *points.get([0, pt]).unwrap();
460            let y = *points.get([1, pt]).unwrap();
461            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0 - x - y);
462            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x);
463            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), y);
464        }
465        check_dofs(e);
466    }
467
468    #[test]
469    fn test_lagrange_0_quadrilateral() {
470        let e = create::<f64, f64>(
471            ReferenceCellType::Quadrilateral,
472            0,
473            Continuity::Discontinuous,
474            LagrangeVariant::Equispaced,
475        );
476        assert_eq!(e.value_size(), 1);
477        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
478        let mut points = rlst_dynamic_array!(f64, [2, 6]);
479        *points.get_mut([0, 0]).unwrap() = 0.0;
480        *points.get_mut([1, 0]).unwrap() = 0.0;
481        *points.get_mut([0, 1]).unwrap() = 1.0;
482        *points.get_mut([1, 1]).unwrap() = 0.0;
483        *points.get_mut([0, 2]).unwrap() = 0.0;
484        *points.get_mut([1, 2]).unwrap() = 1.0;
485        *points.get_mut([0, 3]).unwrap() = 0.5;
486        *points.get_mut([1, 3]).unwrap() = 0.0;
487        *points.get_mut([0, 4]).unwrap() = 0.0;
488        *points.get_mut([1, 4]).unwrap() = 0.5;
489        *points.get_mut([0, 5]).unwrap() = 0.5;
490        *points.get_mut([1, 5]).unwrap() = 0.5;
491        e.tabulate(&points, 0, &mut data);
492
493        for pt in 0..6 {
494            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
495        }
496        check_dofs(e);
497    }
498
499    #[test]
500    fn test_lagrange_1_quadrilateral() {
501        let e = create::<f64, f64>(
502            ReferenceCellType::Quadrilateral,
503            1,
504            Continuity::Standard,
505            LagrangeVariant::Equispaced,
506        );
507        assert_eq!(e.value_size(), 1);
508        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
509        let mut points = rlst_dynamic_array!(f64, [2, 6]);
510        *points.get_mut([0, 0]).unwrap() = 0.0;
511        *points.get_mut([1, 0]).unwrap() = 0.0;
512        *points.get_mut([0, 1]).unwrap() = 1.0;
513        *points.get_mut([1, 1]).unwrap() = 0.0;
514        *points.get_mut([0, 2]).unwrap() = 0.0;
515        *points.get_mut([1, 2]).unwrap() = 1.0;
516        *points.get_mut([0, 3]).unwrap() = 1.0;
517        *points.get_mut([1, 3]).unwrap() = 1.0;
518        *points.get_mut([0, 4]).unwrap() = 0.25;
519        *points.get_mut([1, 4]).unwrap() = 0.5;
520        *points.get_mut([0, 5]).unwrap() = 0.3;
521        *points.get_mut([1, 5]).unwrap() = 0.2;
522
523        e.tabulate(&points, 0, &mut data);
524
525        for pt in 0..6 {
526            let x = *points.get([0, pt]).unwrap();
527            let y = *points.get([1, pt]).unwrap();
528            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), (1.0 - x) * (1.0 - y));
529            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x * (1.0 - y));
530            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), (1.0 - x) * y);
531            assert_relative_eq!(*data.get([0, pt, 3, 0]).unwrap(), x * y);
532        }
533        check_dofs(e);
534    }
535
536    #[test]
537    fn test_lagrange_2_quadrilateral() {
538        let e = create::<f64, f64>(
539            ReferenceCellType::Quadrilateral,
540            2,
541            Continuity::Standard,
542            LagrangeVariant::Equispaced,
543        );
544        assert_eq!(e.value_size(), 1);
545        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
546        let mut points = rlst_dynamic_array!(f64, [2, 6]);
547        *points.get_mut([0, 0]).unwrap() = 0.0;
548        *points.get_mut([1, 0]).unwrap() = 0.0;
549        *points.get_mut([0, 1]).unwrap() = 1.0;
550        *points.get_mut([1, 1]).unwrap() = 0.0;
551        *points.get_mut([0, 2]).unwrap() = 0.0;
552        *points.get_mut([1, 2]).unwrap() = 1.0;
553        *points.get_mut([0, 3]).unwrap() = 1.0;
554        *points.get_mut([1, 3]).unwrap() = 1.0;
555        *points.get_mut([0, 4]).unwrap() = 0.25;
556        *points.get_mut([1, 4]).unwrap() = 0.5;
557        *points.get_mut([0, 5]).unwrap() = 0.3;
558        *points.get_mut([1, 5]).unwrap() = 0.2;
559
560        e.tabulate(&points, 0, &mut data);
561
562        for pt in 0..6 {
563            let x = *points.get([0, pt]).unwrap();
564            let y = *points.get([1, pt]).unwrap();
565            assert_relative_eq!(
566                *data.get([0, pt, 0, 0]).unwrap(),
567                (1.0 - x) * (1.0 - 2.0 * x) * (1.0 - y) * (1.0 - 2.0 * y),
568                epsilon = 1e-14
569            );
570            assert_relative_eq!(
571                *data.get([0, pt, 1, 0]).unwrap(),
572                x * (2.0 * x - 1.0) * (1.0 - y) * (1.0 - 2.0 * y),
573                epsilon = 1e-14
574            );
575            assert_relative_eq!(
576                *data.get([0, pt, 2, 0]).unwrap(),
577                (1.0 - x) * (1.0 - 2.0 * x) * y * (2.0 * y - 1.0),
578                epsilon = 1e-14
579            );
580            assert_relative_eq!(
581                *data.get([0, pt, 3, 0]).unwrap(),
582                x * (2.0 * x - 1.0) * y * (2.0 * y - 1.0),
583                epsilon = 1e-14
584            );
585            assert_relative_eq!(
586                *data.get([0, pt, 4, 0]).unwrap(),
587                4.0 * x * (1.0 - x) * (1.0 - y) * (1.0 - 2.0 * y),
588                epsilon = 1e-14
589            );
590            assert_relative_eq!(
591                *data.get([0, pt, 5, 0]).unwrap(),
592                (1.0 - x) * (1.0 - 2.0 * x) * 4.0 * y * (1.0 - y),
593                epsilon = 1e-14
594            );
595            assert_relative_eq!(
596                *data.get([0, pt, 6, 0]).unwrap(),
597                x * (2.0 * x - 1.0) * 4.0 * y * (1.0 - y),
598                epsilon = 1e-14
599            );
600            assert_relative_eq!(
601                *data.get([0, pt, 7, 0]).unwrap(),
602                4.0 * x * (1.0 - x) * y * (2.0 * y - 1.0),
603                epsilon = 1e-14
604            );
605            assert_relative_eq!(
606                *data.get([0, pt, 8, 0]).unwrap(),
607                4.0 * x * (1.0 - x) * 4.0 * y * (1.0 - y),
608                epsilon = 1e-14
609            );
610        }
611        check_dofs(e);
612    }
613
614    #[test]
615    fn test_lagrange_0_tetrahedron() {
616        let e = create::<f64, f64>(
617            ReferenceCellType::Tetrahedron,
618            0,
619            Continuity::Discontinuous,
620            LagrangeVariant::Equispaced,
621        );
622        assert_eq!(e.value_size(), 1);
623        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
624        let mut points = rlst_dynamic_array!(f64, [3, 6]);
625        *points.get_mut([0, 0]).unwrap() = 0.0;
626        *points.get_mut([1, 0]).unwrap() = 0.0;
627        *points.get_mut([2, 0]).unwrap() = 0.0;
628        *points.get_mut([0, 1]).unwrap() = 1.0;
629        *points.get_mut([1, 1]).unwrap() = 0.0;
630        *points.get_mut([2, 1]).unwrap() = 0.0;
631        *points.get_mut([0, 2]).unwrap() = 0.0;
632        *points.get_mut([1, 2]).unwrap() = 1.0;
633        *points.get_mut([2, 2]).unwrap() = 0.0;
634        *points.get_mut([0, 3]).unwrap() = 0.5;
635        *points.get_mut([1, 3]).unwrap() = 0.0;
636        *points.get_mut([2, 3]).unwrap() = 0.5;
637        *points.get_mut([0, 4]).unwrap() = 0.0;
638        *points.get_mut([1, 4]).unwrap() = 0.5;
639        *points.get_mut([2, 4]).unwrap() = 0.5;
640        *points.get_mut([0, 5]).unwrap() = 0.5;
641        *points.get_mut([1, 5]).unwrap() = 0.2;
642        *points.get_mut([2, 5]).unwrap() = 0.3;
643        e.tabulate(&points, 0, &mut data);
644
645        for pt in 0..6 {
646            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
647        }
648        check_dofs(e);
649    }
650
651    #[test]
652    fn test_lagrange_1_tetrahedron() {
653        let e = create::<f64, f64>(
654            ReferenceCellType::Tetrahedron,
655            1,
656            Continuity::Standard,
657            LagrangeVariant::Equispaced,
658        );
659        assert_eq!(e.value_size(), 1);
660        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
661        let mut points = rlst_dynamic_array!(f64, [3, 6]);
662        *points.get_mut([0, 0]).unwrap() = 0.0;
663        *points.get_mut([1, 0]).unwrap() = 0.0;
664        *points.get_mut([2, 0]).unwrap() = 0.0;
665        *points.get_mut([0, 1]).unwrap() = 1.0;
666        *points.get_mut([1, 1]).unwrap() = 0.0;
667        *points.get_mut([2, 1]).unwrap() = 0.0;
668        *points.get_mut([0, 2]).unwrap() = 0.0;
669        *points.get_mut([1, 2]).unwrap() = 0.8;
670        *points.get_mut([2, 2]).unwrap() = 0.2;
671        *points.get_mut([0, 3]).unwrap() = 0.0;
672        *points.get_mut([1, 3]).unwrap() = 0.0;
673        *points.get_mut([2, 3]).unwrap() = 0.8;
674        *points.get_mut([0, 4]).unwrap() = 0.25;
675        *points.get_mut([1, 4]).unwrap() = 0.5;
676        *points.get_mut([2, 4]).unwrap() = 0.1;
677        *points.get_mut([0, 5]).unwrap() = 0.2;
678        *points.get_mut([1, 5]).unwrap() = 0.1;
679        *points.get_mut([2, 5]).unwrap() = 0.15;
680
681        e.tabulate(&points, 0, &mut data);
682
683        for pt in 0..6 {
684            let x = *points.get([0, pt]).unwrap();
685            let y = *points.get([1, pt]).unwrap();
686            let z = *points.get([2, pt]).unwrap();
687            assert_relative_eq!(
688                *data.get([0, pt, 0, 0]).unwrap(),
689                1.0 - x - y - z,
690                epsilon = 1e-14
691            );
692            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x, epsilon = 1e-14);
693            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), y, epsilon = 1e-14);
694            assert_relative_eq!(*data.get([0, pt, 3, 0]).unwrap(), z, epsilon = 1e-14);
695        }
696        check_dofs(e);
697    }
698
699    #[test]
700    fn test_lagrange_0_hexahedron() {
701        let e = create::<f64, f64>(
702            ReferenceCellType::Hexahedron,
703            0,
704            Continuity::Discontinuous,
705            LagrangeVariant::Equispaced,
706        );
707        assert_eq!(e.value_size(), 1);
708        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
709        let mut points = rlst_dynamic_array!(f64, [3, 6]);
710        *points.get_mut([0, 0]).unwrap() = 0.0;
711        *points.get_mut([1, 0]).unwrap() = 0.0;
712        *points.get_mut([2, 0]).unwrap() = 0.0;
713        *points.get_mut([0, 1]).unwrap() = 1.0;
714        *points.get_mut([1, 1]).unwrap() = 0.0;
715        *points.get_mut([2, 1]).unwrap() = 0.0;
716        *points.get_mut([0, 2]).unwrap() = 0.0;
717        *points.get_mut([1, 2]).unwrap() = 1.0;
718        *points.get_mut([2, 2]).unwrap() = 0.0;
719        *points.get_mut([0, 3]).unwrap() = 0.5;
720        *points.get_mut([1, 3]).unwrap() = 0.0;
721        *points.get_mut([2, 3]).unwrap() = 0.5;
722        *points.get_mut([0, 4]).unwrap() = 0.0;
723        *points.get_mut([1, 4]).unwrap() = 0.5;
724        *points.get_mut([2, 4]).unwrap() = 0.5;
725        *points.get_mut([0, 5]).unwrap() = 0.5;
726        *points.get_mut([1, 5]).unwrap() = 0.5;
727        *points.get_mut([2, 5]).unwrap() = 0.5;
728        e.tabulate(&points, 0, &mut data);
729
730        for pt in 0..6 {
731            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
732        }
733        check_dofs(e);
734    }
735
736    #[test]
737    fn test_lagrange_1_hexahedron() {
738        let e = create::<f64, f64>(
739            ReferenceCellType::Hexahedron,
740            1,
741            Continuity::Standard,
742            LagrangeVariant::Equispaced,
743        );
744        assert_eq!(e.value_size(), 1);
745        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
746        let mut points = rlst_dynamic_array!(f64, [3, 6]);
747        *points.get_mut([0, 0]).unwrap() = 0.0;
748        *points.get_mut([1, 0]).unwrap() = 0.0;
749        *points.get_mut([2, 0]).unwrap() = 0.0;
750        *points.get_mut([0, 1]).unwrap() = 1.0;
751        *points.get_mut([1, 1]).unwrap() = 0.0;
752        *points.get_mut([2, 1]).unwrap() = 0.0;
753        *points.get_mut([0, 2]).unwrap() = 0.0;
754        *points.get_mut([1, 2]).unwrap() = 1.0;
755        *points.get_mut([2, 2]).unwrap() = 1.0;
756        *points.get_mut([0, 3]).unwrap() = 1.0;
757        *points.get_mut([1, 3]).unwrap() = 1.0;
758        *points.get_mut([2, 3]).unwrap() = 1.0;
759        *points.get_mut([0, 4]).unwrap() = 0.25;
760        *points.get_mut([1, 4]).unwrap() = 0.5;
761        *points.get_mut([2, 4]).unwrap() = 0.1;
762        *points.get_mut([0, 5]).unwrap() = 0.3;
763        *points.get_mut([1, 5]).unwrap() = 0.2;
764        *points.get_mut([2, 5]).unwrap() = 0.4;
765
766        e.tabulate(&points, 0, &mut data);
767
768        for pt in 0..6 {
769            let x = *points.get([0, pt]).unwrap();
770            let y = *points.get([1, pt]).unwrap();
771            let z = *points.get([2, pt]).unwrap();
772            assert_relative_eq!(
773                *data.get([0, pt, 0, 0]).unwrap(),
774                (1.0 - x) * (1.0 - y) * (1.0 - z),
775                epsilon = 1e-14
776            );
777            assert_relative_eq!(
778                *data.get([0, pt, 1, 0]).unwrap(),
779                x * (1.0 - y) * (1.0 - z),
780                epsilon = 1e-14
781            );
782            assert_relative_eq!(
783                *data.get([0, pt, 2, 0]).unwrap(),
784                (1.0 - x) * y * (1.0 - z),
785                epsilon = 1e-14
786            );
787            assert_relative_eq!(
788                *data.get([0, pt, 3, 0]).unwrap(),
789                x * y * (1.0 - z),
790                epsilon = 1e-14
791            );
792            assert_relative_eq!(
793                *data.get([0, pt, 4, 0]).unwrap(),
794                (1.0 - x) * (1.0 - y) * z,
795                epsilon = 1e-14
796            );
797            assert_relative_eq!(
798                *data.get([0, pt, 5, 0]).unwrap(),
799                x * (1.0 - y) * z,
800                epsilon = 1e-14
801            );
802            assert_relative_eq!(
803                *data.get([0, pt, 6, 0]).unwrap(),
804                (1.0 - x) * y * z,
805                epsilon = 1e-14
806            );
807            assert_relative_eq!(
808                *data.get([0, pt, 7, 0]).unwrap(),
809                x * y * z,
810                epsilon = 1e-14
811            );
812        }
813        check_dofs(e);
814    }
815
816    #[test]
817    fn test_lagrange_higher_degree_triangle() {
818        create::<f64, f64>(
819            ReferenceCellType::Triangle,
820            2,
821            Continuity::Standard,
822            LagrangeVariant::Equispaced,
823        );
824        create::<f64, f64>(
825            ReferenceCellType::Triangle,
826            3,
827            Continuity::Standard,
828            LagrangeVariant::Equispaced,
829        );
830        create::<f64, f64>(
831            ReferenceCellType::Triangle,
832            4,
833            Continuity::Standard,
834            LagrangeVariant::Equispaced,
835        );
836        create::<f64, f64>(
837            ReferenceCellType::Triangle,
838            5,
839            Continuity::Standard,
840            LagrangeVariant::Equispaced,
841        );
842
843        create::<f64, f64>(
844            ReferenceCellType::Triangle,
845            2,
846            Continuity::Discontinuous,
847            LagrangeVariant::Equispaced,
848        );
849        create::<f64, f64>(
850            ReferenceCellType::Triangle,
851            3,
852            Continuity::Discontinuous,
853            LagrangeVariant::Equispaced,
854        );
855        create::<f64, f64>(
856            ReferenceCellType::Triangle,
857            4,
858            Continuity::Discontinuous,
859            LagrangeVariant::Equispaced,
860        );
861        create::<f64, f64>(
862            ReferenceCellType::Triangle,
863            5,
864            Continuity::Discontinuous,
865            LagrangeVariant::Equispaced,
866        );
867    }
868
869    #[test]
870    fn test_lagrange_higher_degree_interval() {
871        create::<f64, f64>(
872            ReferenceCellType::Interval,
873            2,
874            Continuity::Standard,
875            LagrangeVariant::Equispaced,
876        );
877        create::<f64, f64>(
878            ReferenceCellType::Interval,
879            3,
880            Continuity::Standard,
881            LagrangeVariant::Equispaced,
882        );
883        create::<f64, f64>(
884            ReferenceCellType::Interval,
885            4,
886            Continuity::Standard,
887            LagrangeVariant::Equispaced,
888        );
889        create::<f64, f64>(
890            ReferenceCellType::Interval,
891            5,
892            Continuity::Standard,
893            LagrangeVariant::Equispaced,
894        );
895
896        create::<f64, f64>(
897            ReferenceCellType::Interval,
898            2,
899            Continuity::Discontinuous,
900            LagrangeVariant::Equispaced,
901        );
902        create::<f64, f64>(
903            ReferenceCellType::Interval,
904            3,
905            Continuity::Discontinuous,
906            LagrangeVariant::Equispaced,
907        );
908        create::<f64, f64>(
909            ReferenceCellType::Interval,
910            4,
911            Continuity::Discontinuous,
912            LagrangeVariant::Equispaced,
913        );
914        create::<f64, f64>(
915            ReferenceCellType::Interval,
916            5,
917            Continuity::Discontinuous,
918            LagrangeVariant::Equispaced,
919        );
920    }
921
922    #[test]
923    fn test_lagrange_higher_degree_quadrilateral() {
924        create::<f64, f64>(
925            ReferenceCellType::Quadrilateral,
926            2,
927            Continuity::Standard,
928            LagrangeVariant::Equispaced,
929        );
930        create::<f64, f64>(
931            ReferenceCellType::Quadrilateral,
932            3,
933            Continuity::Standard,
934            LagrangeVariant::Equispaced,
935        );
936        create::<f64, f64>(
937            ReferenceCellType::Quadrilateral,
938            4,
939            Continuity::Standard,
940            LagrangeVariant::Equispaced,
941        );
942        create::<f64, f64>(
943            ReferenceCellType::Quadrilateral,
944            5,
945            Continuity::Standard,
946            LagrangeVariant::Equispaced,
947        );
948
949        create::<f64, f64>(
950            ReferenceCellType::Quadrilateral,
951            2,
952            Continuity::Discontinuous,
953            LagrangeVariant::Equispaced,
954        );
955        create::<f64, f64>(
956            ReferenceCellType::Quadrilateral,
957            3,
958            Continuity::Discontinuous,
959            LagrangeVariant::Equispaced,
960        );
961        create::<f64, f64>(
962            ReferenceCellType::Quadrilateral,
963            4,
964            Continuity::Discontinuous,
965            LagrangeVariant::Equispaced,
966        );
967        create::<f64, f64>(
968            ReferenceCellType::Quadrilateral,
969            5,
970            Continuity::Discontinuous,
971            LagrangeVariant::Equispaced,
972        );
973    }
974
975    #[test]
976    fn test_lagrange_interval_equispaced() {
977        let e = create::<f64, f64>(
978            ReferenceCellType::Interval,
979            4,
980            Continuity::Discontinuous,
981            LagrangeVariant::Equispaced,
982        );
983
984        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 2));
985        let mut points = rlst_dynamic_array!(f64, [1, 2]);
986        *points.get_mut([0, 0]).unwrap() = 0.25;
987        *points.get_mut([0, 1]).unwrap() = 0.75;
988
989        e.tabulate(&points, 0, &mut data);
990
991        for pt in 0..data.shape()[1] {
992            for basis in 0..data.shape()[1] {
993                assert!(
994                    data.get([0, pt, basis, 0]).unwrap().abs() < 1e-10
995                        || (1.0 - data.get([0, pt, basis, 0]).unwrap()).abs() < 1e-14
996                );
997            }
998        }
999    }
1000
1001    #[test]
1002    fn test_lagrange_interval_gll() {
1003        let e = create::<f64, f64>(
1004            ReferenceCellType::Interval,
1005            4,
1006            Continuity::Discontinuous,
1007            LagrangeVariant::GLL,
1008        );
1009
1010        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 2));
1011        let mut points = rlst_dynamic_array!(f64, [1, 2]);
1012        *points.get_mut([0, 0]).unwrap() = 0.25;
1013        *points.get_mut([0, 1]).unwrap() = 0.75;
1014
1015        e.tabulate(&points, 0, &mut data);
1016
1017        for pt in 0..data.shape()[1] {
1018            for basis in 0..data.shape()[1] {
1019                assert!(data.get([0, pt, basis, 0]).unwrap().abs() > 1e-10);
1020            }
1021        }
1022    }
1023
1024    #[test]
1025    fn test_lagrange_quadrilateral_equispaced() {
1026        let e = create::<f64, f64>(
1027            ReferenceCellType::Quadrilateral,
1028            4,
1029            Continuity::Discontinuous,
1030            LagrangeVariant::Equispaced,
1031        );
1032
1033        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
1034        let mut points = rlst_dynamic_array!(f64, [2, 4]);
1035        *points.get_mut([0, 0]).unwrap() = 0.25;
1036        *points.get_mut([1, 0]).unwrap() = 0.25;
1037        *points.get_mut([0, 1]).unwrap() = 0.25;
1038        *points.get_mut([1, 1]).unwrap() = 0.75;
1039        *points.get_mut([0, 2]).unwrap() = 0.75;
1040        *points.get_mut([1, 2]).unwrap() = 0.25;
1041        *points.get_mut([0, 3]).unwrap() = 0.75;
1042        *points.get_mut([1, 3]).unwrap() = 0.75;
1043
1044        e.tabulate(&points, 0, &mut data);
1045
1046        for pt in 0..data.shape()[1] {
1047            for basis in 0..data.shape()[1] {
1048                assert!(
1049                    data.get([0, pt, basis, 0]).unwrap().abs() < 1e-10
1050                        || (1.0 - data.get([0, pt, basis, 0]).unwrap()).abs() < 1e-14
1051                );
1052            }
1053        }
1054    }
1055
1056    #[test]
1057    fn test_lagrange_quadrilateral_gll() {
1058        let e = create::<f64, f64>(
1059            ReferenceCellType::Quadrilateral,
1060            4,
1061            Continuity::Discontinuous,
1062            LagrangeVariant::GLL,
1063        );
1064
1065        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
1066        let mut points = rlst_dynamic_array!(f64, [2, 4]);
1067        *points.get_mut([0, 0]).unwrap() = 0.25;
1068        *points.get_mut([1, 0]).unwrap() = 0.25;
1069        *points.get_mut([0, 1]).unwrap() = 0.25;
1070        *points.get_mut([1, 1]).unwrap() = 0.75;
1071        *points.get_mut([0, 2]).unwrap() = 0.75;
1072        *points.get_mut([1, 2]).unwrap() = 0.25;
1073        *points.get_mut([0, 3]).unwrap() = 0.75;
1074        *points.get_mut([1, 3]).unwrap() = 0.75;
1075
1076        e.tabulate(&points, 0, &mut data);
1077
1078        for pt in 0..data.shape()[1] {
1079            for basis in 0..data.shape()[1] {
1080                assert!(data.get([0, pt, basis, 0]).unwrap().abs() > 1e-10);
1081            }
1082        }
1083    }
1084
1085    fn ordered(x: &[f64], y: &[f64]) -> bool {
1086        assert_eq!(x.len(), y.len());
1087        if x.is_empty() {
1088            false
1089        } else if (x[x.len() - 1] - y[y.len() - 1]).abs() < 1e-8 {
1090            ordered(&x[..x.len() - 1], &y[..y.len() - 1])
1091        } else {
1092            x[x.len() - 1] < y[y.len() - 1]
1093        }
1094    }
1095
1096    macro_rules! test_point_ordering {
1097        ($cell:ident, $max_degree:expr) => {
1098            paste! {
1099                #[test]
1100                fn [<test_point_ordering_ $cell:lower>]() {
1101                    for degree in 3..=[<$max_degree>] {
1102                        let e = create::<f64, f64>(
1103                            ReferenceCellType::[<$cell>],
1104                            degree,
1105                            Continuity::Standard,
1106                            LagrangeVariant::Equispaced,
1107                        );
1108                        for pt_sets in &e.interpolation_points {
1109                            for pts in pt_sets {
1110                                if pts.shape()[1] > 0 {
1111                                    for i in 0..pts.shape()[1] - 1 {
1112                                        assert!(ordered(&pts.r().col(i).data().unwrap(), &pts.r().col(i + 1).data().unwrap()));
1113                                    }
1114                                }
1115                            }
1116                        }
1117                    }
1118                }
1119            }
1120        };
1121    }
1122
1123    test_point_ordering!(Triangle, 7);
1124    test_point_ordering!(Quadrilateral, 7);
1125    test_point_ordering!(Tetrahedron, 6);
1126    test_point_ordering!(Hexahedron, 4);
1127}
ndelement/ciarlet/lagrange.rs

ndelement/ciarlet/
lagrange.rs
