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 Variant {
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: Variant,
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            Variant::Equispaced => (1..degree)
82                .map(|i| num::cast::<_, TGeo>(i).unwrap() / num::cast::<_, TGeo>(degree).unwrap())
83                .collect::<Vec<_>>(),
84            Variant::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 != Variant::Equispaced {
145                        unimplemented!();
146                    }
147                    let mut n = 0;
148                    for i0 in 1..degree {
149                        for i1 in 1..degree - i0 {
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 p0 in &pts1d {
167                        for p1 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 != Variant::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 i0 in 1..degree {
224                        for i1 in 1..degree - i0 {
225                            for i2 in 1..degree - i0 - 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 p0 in &pts1d {
255                        for p1 in &pts1d {
256                            for p2 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: Variant,
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: Variant) -> 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 rlst::DynArray;
339
340    #[test]
341    fn test_lagrange_1() {
342        let e = create::<f64, f64>(
343            ReferenceCellType::Triangle,
344            1,
345            Continuity::Standard,
346            Variant::Equispaced,
347        );
348        assert_eq!(e.value_size(), 1);
349    }
350
351    #[test]
352    fn test_lagrange_0_interval() {
353        let e = create::<f64, f64>(
354            ReferenceCellType::Interval,
355            0,
356            Continuity::Discontinuous,
357            Variant::Equispaced,
358        );
359        assert_eq!(e.value_size(), 1);
360        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
361        let mut points = rlst_dynamic_array!(f64, [1, 4]);
362        *points.get_mut([0, 0]).unwrap() = 0.0;
363        *points.get_mut([0, 1]).unwrap() = 0.2;
364        *points.get_mut([0, 2]).unwrap() = 0.4;
365        *points.get_mut([0, 3]).unwrap() = 1.0;
366        e.tabulate(&points, 0, &mut data);
367
368        for pt in 0..4 {
369            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
370        }
371        check_dofs(e);
372    }
373
374    #[test]
375    fn test_lagrange_1_interval() {
376        let e = create::<f64, f64>(
377            ReferenceCellType::Interval,
378            1,
379            Continuity::Standard,
380            Variant::Equispaced,
381        );
382        assert_eq!(e.value_size(), 1);
383        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
384        let mut points = rlst_dynamic_array!(f64, [1, 4]);
385        *points.get_mut([0, 0]).unwrap() = 0.0;
386        *points.get_mut([0, 1]).unwrap() = 0.2;
387        *points.get_mut([0, 2]).unwrap() = 0.4;
388        *points.get_mut([0, 3]).unwrap() = 1.0;
389        e.tabulate(&points, 0, &mut data);
390
391        for pt in 0..4 {
392            let x = *points.get([0, pt]).unwrap();
393            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0 - x);
394            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x);
395        }
396        check_dofs(e);
397    }
398
399    #[test]
400    fn test_lagrange_0_triangle() {
401        let e = create::<f64, f64>(
402            ReferenceCellType::Triangle,
403            0,
404            Continuity::Discontinuous,
405            Variant::Equispaced,
406        );
407        assert_eq!(e.value_size(), 1);
408        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
409
410        let mut points = rlst_dynamic_array!(f64, [2, 6]);
411        *points.get_mut([0, 0]).unwrap() = 0.0;
412        *points.get_mut([1, 0]).unwrap() = 0.0;
413        *points.get_mut([0, 1]).unwrap() = 1.0;
414        *points.get_mut([1, 1]).unwrap() = 0.0;
415        *points.get_mut([0, 2]).unwrap() = 0.0;
416        *points.get_mut([1, 2]).unwrap() = 1.0;
417        *points.get_mut([0, 3]).unwrap() = 0.5;
418        *points.get_mut([1, 3]).unwrap() = 0.0;
419        *points.get_mut([0, 4]).unwrap() = 0.0;
420        *points.get_mut([1, 4]).unwrap() = 0.5;
421        *points.get_mut([0, 5]).unwrap() = 0.5;
422        *points.get_mut([1, 5]).unwrap() = 0.5;
423
424        e.tabulate(&points, 0, &mut data);
425
426        for pt in 0..6 {
427            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
428        }
429        check_dofs(e);
430    }
431
432    #[test]
433    fn test_lagrange_1_triangle() {
434        let e = create::<f64, f64>(
435            ReferenceCellType::Triangle,
436            1,
437            Continuity::Standard,
438            Variant::Equispaced,
439        );
440        assert_eq!(e.value_size(), 1);
441        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
442        let mut points = rlst_dynamic_array!(f64, [2, 6]);
443        *points.get_mut([0, 0]).unwrap() = 0.0;
444        *points.get_mut([1, 0]).unwrap() = 0.0;
445        *points.get_mut([0, 1]).unwrap() = 1.0;
446        *points.get_mut([1, 1]).unwrap() = 0.0;
447        *points.get_mut([0, 2]).unwrap() = 0.0;
448        *points.get_mut([1, 2]).unwrap() = 1.0;
449        *points.get_mut([0, 3]).unwrap() = 0.5;
450        *points.get_mut([1, 3]).unwrap() = 0.0;
451        *points.get_mut([0, 4]).unwrap() = 0.0;
452        *points.get_mut([1, 4]).unwrap() = 0.5;
453        *points.get_mut([0, 5]).unwrap() = 0.5;
454        *points.get_mut([1, 5]).unwrap() = 0.5;
455        e.tabulate(&points, 0, &mut data);
456
457        for pt in 0..6 {
458            let x = *points.get([0, pt]).unwrap();
459            let y = *points.get([1, pt]).unwrap();
460            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0 - x - y);
461            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x);
462            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), y);
463        }
464        check_dofs(e);
465    }
466
467    #[test]
468    fn test_lagrange_0_quadrilateral() {
469        let e = create::<f64, f64>(
470            ReferenceCellType::Quadrilateral,
471            0,
472            Continuity::Discontinuous,
473            Variant::Equispaced,
474        );
475        assert_eq!(e.value_size(), 1);
476        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
477        let mut points = rlst_dynamic_array!(f64, [2, 6]);
478        *points.get_mut([0, 0]).unwrap() = 0.0;
479        *points.get_mut([1, 0]).unwrap() = 0.0;
480        *points.get_mut([0, 1]).unwrap() = 1.0;
481        *points.get_mut([1, 1]).unwrap() = 0.0;
482        *points.get_mut([0, 2]).unwrap() = 0.0;
483        *points.get_mut([1, 2]).unwrap() = 1.0;
484        *points.get_mut([0, 3]).unwrap() = 0.5;
485        *points.get_mut([1, 3]).unwrap() = 0.0;
486        *points.get_mut([0, 4]).unwrap() = 0.0;
487        *points.get_mut([1, 4]).unwrap() = 0.5;
488        *points.get_mut([0, 5]).unwrap() = 0.5;
489        *points.get_mut([1, 5]).unwrap() = 0.5;
490        e.tabulate(&points, 0, &mut data);
491
492        for pt in 0..6 {
493            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
494        }
495        check_dofs(e);
496    }
497
498    #[test]
499    fn test_lagrange_1_quadrilateral() {
500        let e = create::<f64, f64>(
501            ReferenceCellType::Quadrilateral,
502            1,
503            Continuity::Standard,
504            Variant::Equispaced,
505        );
506        assert_eq!(e.value_size(), 1);
507        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
508        let mut points = rlst_dynamic_array!(f64, [2, 6]);
509        *points.get_mut([0, 0]).unwrap() = 0.0;
510        *points.get_mut([1, 0]).unwrap() = 0.0;
511        *points.get_mut([0, 1]).unwrap() = 1.0;
512        *points.get_mut([1, 1]).unwrap() = 0.0;
513        *points.get_mut([0, 2]).unwrap() = 0.0;
514        *points.get_mut([1, 2]).unwrap() = 1.0;
515        *points.get_mut([0, 3]).unwrap() = 1.0;
516        *points.get_mut([1, 3]).unwrap() = 1.0;
517        *points.get_mut([0, 4]).unwrap() = 0.25;
518        *points.get_mut([1, 4]).unwrap() = 0.5;
519        *points.get_mut([0, 5]).unwrap() = 0.3;
520        *points.get_mut([1, 5]).unwrap() = 0.2;
521
522        e.tabulate(&points, 0, &mut data);
523
524        for pt in 0..6 {
525            let x = *points.get([0, pt]).unwrap();
526            let y = *points.get([1, pt]).unwrap();
527            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), (1.0 - x) * (1.0 - y));
528            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x * (1.0 - y));
529            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), (1.0 - x) * y);
530            assert_relative_eq!(*data.get([0, pt, 3, 0]).unwrap(), x * y);
531        }
532        check_dofs(e);
533    }
534
535    #[test]
536    fn test_lagrange_2_quadrilateral() {
537        let e = create::<f64, f64>(
538            ReferenceCellType::Quadrilateral,
539            2,
540            Continuity::Standard,
541            Variant::Equispaced,
542        );
543        assert_eq!(e.value_size(), 1);
544        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
545        let mut points = rlst_dynamic_array!(f64, [2, 6]);
546        *points.get_mut([0, 0]).unwrap() = 0.0;
547        *points.get_mut([1, 0]).unwrap() = 0.0;
548        *points.get_mut([0, 1]).unwrap() = 1.0;
549        *points.get_mut([1, 1]).unwrap() = 0.0;
550        *points.get_mut([0, 2]).unwrap() = 0.0;
551        *points.get_mut([1, 2]).unwrap() = 1.0;
552        *points.get_mut([0, 3]).unwrap() = 1.0;
553        *points.get_mut([1, 3]).unwrap() = 1.0;
554        *points.get_mut([0, 4]).unwrap() = 0.25;
555        *points.get_mut([1, 4]).unwrap() = 0.5;
556        *points.get_mut([0, 5]).unwrap() = 0.3;
557        *points.get_mut([1, 5]).unwrap() = 0.2;
558
559        e.tabulate(&points, 0, &mut data);
560
561        for pt in 0..6 {
562            let x = *points.get([0, pt]).unwrap();
563            let y = *points.get([1, pt]).unwrap();
564            assert_relative_eq!(
565                *data.get([0, pt, 0, 0]).unwrap(),
566                (1.0 - x) * (1.0 - 2.0 * x) * (1.0 - y) * (1.0 - 2.0 * y),
567                epsilon = 1e-14
568            );
569            assert_relative_eq!(
570                *data.get([0, pt, 1, 0]).unwrap(),
571                x * (2.0 * x - 1.0) * (1.0 - y) * (1.0 - 2.0 * y),
572                epsilon = 1e-14
573            );
574            assert_relative_eq!(
575                *data.get([0, pt, 2, 0]).unwrap(),
576                (1.0 - x) * (1.0 - 2.0 * x) * y * (2.0 * y - 1.0),
577                epsilon = 1e-14
578            );
579            assert_relative_eq!(
580                *data.get([0, pt, 3, 0]).unwrap(),
581                x * (2.0 * x - 1.0) * y * (2.0 * y - 1.0),
582                epsilon = 1e-14
583            );
584            assert_relative_eq!(
585                *data.get([0, pt, 4, 0]).unwrap(),
586                4.0 * x * (1.0 - x) * (1.0 - y) * (1.0 - 2.0 * y),
587                epsilon = 1e-14
588            );
589            assert_relative_eq!(
590                *data.get([0, pt, 5, 0]).unwrap(),
591                (1.0 - x) * (1.0 - 2.0 * x) * 4.0 * y * (1.0 - y),
592                epsilon = 1e-14
593            );
594            assert_relative_eq!(
595                *data.get([0, pt, 6, 0]).unwrap(),
596                x * (2.0 * x - 1.0) * 4.0 * y * (1.0 - y),
597                epsilon = 1e-14
598            );
599            assert_relative_eq!(
600                *data.get([0, pt, 7, 0]).unwrap(),
601                4.0 * x * (1.0 - x) * y * (2.0 * y - 1.0),
602                epsilon = 1e-14
603            );
604            assert_relative_eq!(
605                *data.get([0, pt, 8, 0]).unwrap(),
606                4.0 * x * (1.0 - x) * 4.0 * y * (1.0 - y),
607                epsilon = 1e-14
608            );
609        }
610        check_dofs(e);
611    }
612
613    #[test]
614    fn test_lagrange_0_tetrahedron() {
615        let e = create::<f64, f64>(
616            ReferenceCellType::Tetrahedron,
617            0,
618            Continuity::Discontinuous,
619            Variant::Equispaced,
620        );
621        assert_eq!(e.value_size(), 1);
622        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
623        let mut points = rlst_dynamic_array!(f64, [3, 6]);
624        *points.get_mut([0, 0]).unwrap() = 0.0;
625        *points.get_mut([1, 0]).unwrap() = 0.0;
626        *points.get_mut([2, 0]).unwrap() = 0.0;
627        *points.get_mut([0, 1]).unwrap() = 1.0;
628        *points.get_mut([1, 1]).unwrap() = 0.0;
629        *points.get_mut([2, 1]).unwrap() = 0.0;
630        *points.get_mut([0, 2]).unwrap() = 0.0;
631        *points.get_mut([1, 2]).unwrap() = 1.0;
632        *points.get_mut([2, 2]).unwrap() = 0.0;
633        *points.get_mut([0, 3]).unwrap() = 0.5;
634        *points.get_mut([1, 3]).unwrap() = 0.0;
635        *points.get_mut([2, 3]).unwrap() = 0.5;
636        *points.get_mut([0, 4]).unwrap() = 0.0;
637        *points.get_mut([1, 4]).unwrap() = 0.5;
638        *points.get_mut([2, 4]).unwrap() = 0.5;
639        *points.get_mut([0, 5]).unwrap() = 0.5;
640        *points.get_mut([1, 5]).unwrap() = 0.2;
641        *points.get_mut([2, 5]).unwrap() = 0.3;
642        e.tabulate(&points, 0, &mut data);
643
644        for pt in 0..6 {
645            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
646        }
647        check_dofs(e);
648    }
649
650    #[test]
651    fn test_lagrange_1_tetrahedron() {
652        let e = create::<f64, f64>(
653            ReferenceCellType::Tetrahedron,
654            1,
655            Continuity::Standard,
656            Variant::Equispaced,
657        );
658        assert_eq!(e.value_size(), 1);
659        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
660        let mut points = rlst_dynamic_array!(f64, [3, 6]);
661        *points.get_mut([0, 0]).unwrap() = 0.0;
662        *points.get_mut([1, 0]).unwrap() = 0.0;
663        *points.get_mut([2, 0]).unwrap() = 0.0;
664        *points.get_mut([0, 1]).unwrap() = 1.0;
665        *points.get_mut([1, 1]).unwrap() = 0.0;
666        *points.get_mut([2, 1]).unwrap() = 0.0;
667        *points.get_mut([0, 2]).unwrap() = 0.0;
668        *points.get_mut([1, 2]).unwrap() = 0.8;
669        *points.get_mut([2, 2]).unwrap() = 0.2;
670        *points.get_mut([0, 3]).unwrap() = 0.0;
671        *points.get_mut([1, 3]).unwrap() = 0.0;
672        *points.get_mut([2, 3]).unwrap() = 0.8;
673        *points.get_mut([0, 4]).unwrap() = 0.25;
674        *points.get_mut([1, 4]).unwrap() = 0.5;
675        *points.get_mut([2, 4]).unwrap() = 0.1;
676        *points.get_mut([0, 5]).unwrap() = 0.2;
677        *points.get_mut([1, 5]).unwrap() = 0.1;
678        *points.get_mut([2, 5]).unwrap() = 0.15;
679
680        e.tabulate(&points, 0, &mut data);
681
682        for pt in 0..6 {
683            let x = *points.get([0, pt]).unwrap();
684            let y = *points.get([1, pt]).unwrap();
685            let z = *points.get([2, pt]).unwrap();
686            assert_relative_eq!(
687                *data.get([0, pt, 0, 0]).unwrap(),
688                1.0 - x - y - z,
689                epsilon = 1e-14
690            );
691            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x, epsilon = 1e-14);
692            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), y, epsilon = 1e-14);
693            assert_relative_eq!(*data.get([0, pt, 3, 0]).unwrap(), z, epsilon = 1e-14);
694        }
695        check_dofs(e);
696    }
697
698    #[test]
699    fn test_lagrange_0_hexahedron() {
700        let e = create::<f64, f64>(
701            ReferenceCellType::Hexahedron,
702            0,
703            Continuity::Discontinuous,
704            Variant::Equispaced,
705        );
706        assert_eq!(e.value_size(), 1);
707        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
708        let mut points = rlst_dynamic_array!(f64, [3, 6]);
709        *points.get_mut([0, 0]).unwrap() = 0.0;
710        *points.get_mut([1, 0]).unwrap() = 0.0;
711        *points.get_mut([2, 0]).unwrap() = 0.0;
712        *points.get_mut([0, 1]).unwrap() = 1.0;
713        *points.get_mut([1, 1]).unwrap() = 0.0;
714        *points.get_mut([2, 1]).unwrap() = 0.0;
715        *points.get_mut([0, 2]).unwrap() = 0.0;
716        *points.get_mut([1, 2]).unwrap() = 1.0;
717        *points.get_mut([2, 2]).unwrap() = 0.0;
718        *points.get_mut([0, 3]).unwrap() = 0.5;
719        *points.get_mut([1, 3]).unwrap() = 0.0;
720        *points.get_mut([2, 3]).unwrap() = 0.5;
721        *points.get_mut([0, 4]).unwrap() = 0.0;
722        *points.get_mut([1, 4]).unwrap() = 0.5;
723        *points.get_mut([2, 4]).unwrap() = 0.5;
724        *points.get_mut([0, 5]).unwrap() = 0.5;
725        *points.get_mut([1, 5]).unwrap() = 0.5;
726        *points.get_mut([2, 5]).unwrap() = 0.5;
727        e.tabulate(&points, 0, &mut data);
728
729        for pt in 0..6 {
730            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0);
731        }
732        check_dofs(e);
733    }
734
735    #[test]
736    fn test_lagrange_1_hexahedron() {
737        let e = create::<f64, f64>(
738            ReferenceCellType::Hexahedron,
739            1,
740            Continuity::Standard,
741            Variant::Equispaced,
742        );
743        assert_eq!(e.value_size(), 1);
744        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 6));
745        let mut points = rlst_dynamic_array!(f64, [3, 6]);
746        *points.get_mut([0, 0]).unwrap() = 0.0;
747        *points.get_mut([1, 0]).unwrap() = 0.0;
748        *points.get_mut([2, 0]).unwrap() = 0.0;
749        *points.get_mut([0, 1]).unwrap() = 1.0;
750        *points.get_mut([1, 1]).unwrap() = 0.0;
751        *points.get_mut([2, 1]).unwrap() = 0.0;
752        *points.get_mut([0, 2]).unwrap() = 0.0;
753        *points.get_mut([1, 2]).unwrap() = 1.0;
754        *points.get_mut([2, 2]).unwrap() = 1.0;
755        *points.get_mut([0, 3]).unwrap() = 1.0;
756        *points.get_mut([1, 3]).unwrap() = 1.0;
757        *points.get_mut([2, 3]).unwrap() = 1.0;
758        *points.get_mut([0, 4]).unwrap() = 0.25;
759        *points.get_mut([1, 4]).unwrap() = 0.5;
760        *points.get_mut([2, 4]).unwrap() = 0.1;
761        *points.get_mut([0, 5]).unwrap() = 0.3;
762        *points.get_mut([1, 5]).unwrap() = 0.2;
763        *points.get_mut([2, 5]).unwrap() = 0.4;
764
765        e.tabulate(&points, 0, &mut data);
766
767        for pt in 0..6 {
768            let x = *points.get([0, pt]).unwrap();
769            let y = *points.get([1, pt]).unwrap();
770            let z = *points.get([2, pt]).unwrap();
771            assert_relative_eq!(
772                *data.get([0, pt, 0, 0]).unwrap(),
773                (1.0 - x) * (1.0 - y) * (1.0 - z),
774                epsilon = 1e-14
775            );
776            assert_relative_eq!(
777                *data.get([0, pt, 1, 0]).unwrap(),
778                x * (1.0 - y) * (1.0 - z),
779                epsilon = 1e-14
780            );
781            assert_relative_eq!(
782                *data.get([0, pt, 2, 0]).unwrap(),
783                (1.0 - x) * y * (1.0 - z),
784                epsilon = 1e-14
785            );
786            assert_relative_eq!(
787                *data.get([0, pt, 3, 0]).unwrap(),
788                x * y * (1.0 - z),
789                epsilon = 1e-14
790            );
791            assert_relative_eq!(
792                *data.get([0, pt, 4, 0]).unwrap(),
793                (1.0 - x) * (1.0 - y) * z,
794                epsilon = 1e-14
795            );
796            assert_relative_eq!(
797                *data.get([0, pt, 5, 0]).unwrap(),
798                x * (1.0 - y) * z,
799                epsilon = 1e-14
800            );
801            assert_relative_eq!(
802                *data.get([0, pt, 6, 0]).unwrap(),
803                (1.0 - x) * y * z,
804                epsilon = 1e-14
805            );
806            assert_relative_eq!(
807                *data.get([0, pt, 7, 0]).unwrap(),
808                x * y * z,
809                epsilon = 1e-14
810            );
811        }
812        check_dofs(e);
813    }
814
815    #[test]
816    fn test_lagrange_higher_degree_triangle() {
817        create::<f64, f64>(
818            ReferenceCellType::Triangle,
819            2,
820            Continuity::Standard,
821            Variant::Equispaced,
822        );
823        create::<f64, f64>(
824            ReferenceCellType::Triangle,
825            3,
826            Continuity::Standard,
827            Variant::Equispaced,
828        );
829        create::<f64, f64>(
830            ReferenceCellType::Triangle,
831            4,
832            Continuity::Standard,
833            Variant::Equispaced,
834        );
835        create::<f64, f64>(
836            ReferenceCellType::Triangle,
837            5,
838            Continuity::Standard,
839            Variant::Equispaced,
840        );
841
842        create::<f64, f64>(
843            ReferenceCellType::Triangle,
844            2,
845            Continuity::Discontinuous,
846            Variant::Equispaced,
847        );
848        create::<f64, f64>(
849            ReferenceCellType::Triangle,
850            3,
851            Continuity::Discontinuous,
852            Variant::Equispaced,
853        );
854        create::<f64, f64>(
855            ReferenceCellType::Triangle,
856            4,
857            Continuity::Discontinuous,
858            Variant::Equispaced,
859        );
860        create::<f64, f64>(
861            ReferenceCellType::Triangle,
862            5,
863            Continuity::Discontinuous,
864            Variant::Equispaced,
865        );
866    }
867
868    #[test]
869    fn test_lagrange_higher_degree_interval() {
870        create::<f64, f64>(
871            ReferenceCellType::Interval,
872            2,
873            Continuity::Standard,
874            Variant::Equispaced,
875        );
876        create::<f64, f64>(
877            ReferenceCellType::Interval,
878            3,
879            Continuity::Standard,
880            Variant::Equispaced,
881        );
882        create::<f64, f64>(
883            ReferenceCellType::Interval,
884            4,
885            Continuity::Standard,
886            Variant::Equispaced,
887        );
888        create::<f64, f64>(
889            ReferenceCellType::Interval,
890            5,
891            Continuity::Standard,
892            Variant::Equispaced,
893        );
894
895        create::<f64, f64>(
896            ReferenceCellType::Interval,
897            2,
898            Continuity::Discontinuous,
899            Variant::Equispaced,
900        );
901        create::<f64, f64>(
902            ReferenceCellType::Interval,
903            3,
904            Continuity::Discontinuous,
905            Variant::Equispaced,
906        );
907        create::<f64, f64>(
908            ReferenceCellType::Interval,
909            4,
910            Continuity::Discontinuous,
911            Variant::Equispaced,
912        );
913        create::<f64, f64>(
914            ReferenceCellType::Interval,
915            5,
916            Continuity::Discontinuous,
917            Variant::Equispaced,
918        );
919    }
920
921    #[test]
922    fn test_lagrange_higher_degree_quadrilateral() {
923        create::<f64, f64>(
924            ReferenceCellType::Quadrilateral,
925            2,
926            Continuity::Standard,
927            Variant::Equispaced,
928        );
929        create::<f64, f64>(
930            ReferenceCellType::Quadrilateral,
931            3,
932            Continuity::Standard,
933            Variant::Equispaced,
934        );
935        create::<f64, f64>(
936            ReferenceCellType::Quadrilateral,
937            4,
938            Continuity::Standard,
939            Variant::Equispaced,
940        );
941        create::<f64, f64>(
942            ReferenceCellType::Quadrilateral,
943            5,
944            Continuity::Standard,
945            Variant::Equispaced,
946        );
947
948        create::<f64, f64>(
949            ReferenceCellType::Quadrilateral,
950            2,
951            Continuity::Discontinuous,
952            Variant::Equispaced,
953        );
954        create::<f64, f64>(
955            ReferenceCellType::Quadrilateral,
956            3,
957            Continuity::Discontinuous,
958            Variant::Equispaced,
959        );
960        create::<f64, f64>(
961            ReferenceCellType::Quadrilateral,
962            4,
963            Continuity::Discontinuous,
964            Variant::Equispaced,
965        );
966        create::<f64, f64>(
967            ReferenceCellType::Quadrilateral,
968            5,
969            Continuity::Discontinuous,
970            Variant::Equispaced,
971        );
972    }
973
974    #[test]
975    fn test_lagrange_interval_equispaced() {
976        let e = create::<f64, f64>(
977            ReferenceCellType::Interval,
978            4,
979            Continuity::Discontinuous,
980            Variant::Equispaced,
981        );
982
983        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 2));
984        let mut points = rlst_dynamic_array!(f64, [1, 2]);
985        *points.get_mut([0, 0]).unwrap() = 0.25;
986        *points.get_mut([0, 1]).unwrap() = 0.75;
987
988        e.tabulate(&points, 0, &mut data);
989
990        for pt in 0..data.shape()[1] {
991            for basis in 0..data.shape()[1] {
992                assert!(
993                    data.get([0, pt, basis, 0]).unwrap().abs() < 1e-10
994                        || (1.0 - data.get([0, pt, basis, 0]).unwrap()).abs() < 1e-14
995                );
996            }
997        }
998    }
999
1000    #[test]
1001    fn test_lagrange_interval_gll() {
1002        let e = create::<f64, f64>(
1003            ReferenceCellType::Interval,
1004            4,
1005            Continuity::Discontinuous,
1006            Variant::GLL,
1007        );
1008
1009        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 2));
1010        let mut points = rlst_dynamic_array!(f64, [1, 2]);
1011        *points.get_mut([0, 0]).unwrap() = 0.25;
1012        *points.get_mut([0, 1]).unwrap() = 0.75;
1013
1014        e.tabulate(&points, 0, &mut data);
1015
1016        for pt in 0..data.shape()[1] {
1017            for basis in 0..data.shape()[1] {
1018                assert!(data.get([0, pt, basis, 0]).unwrap().abs() > 1e-10);
1019            }
1020        }
1021    }
1022
1023    #[test]
1024    fn test_lagrange_quadrilateral_equispaced() {
1025        let e = create::<f64, f64>(
1026            ReferenceCellType::Quadrilateral,
1027            4,
1028            Continuity::Discontinuous,
1029            Variant::Equispaced,
1030        );
1031
1032        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
1033        let mut points = rlst_dynamic_array!(f64, [2, 4]);
1034        *points.get_mut([0, 0]).unwrap() = 0.25;
1035        *points.get_mut([1, 0]).unwrap() = 0.25;
1036        *points.get_mut([0, 1]).unwrap() = 0.25;
1037        *points.get_mut([1, 1]).unwrap() = 0.75;
1038        *points.get_mut([0, 2]).unwrap() = 0.75;
1039        *points.get_mut([1, 2]).unwrap() = 0.25;
1040        *points.get_mut([0, 3]).unwrap() = 0.75;
1041        *points.get_mut([1, 3]).unwrap() = 0.75;
1042
1043        e.tabulate(&points, 0, &mut data);
1044
1045        for pt in 0..data.shape()[1] {
1046            for basis in 0..data.shape()[1] {
1047                assert!(
1048                    data.get([0, pt, basis, 0]).unwrap().abs() < 1e-10
1049                        || (1.0 - data.get([0, pt, basis, 0]).unwrap()).abs() < 1e-14
1050                );
1051            }
1052        }
1053    }
1054
1055    #[test]
1056    fn test_lagrange_quadrilateral_gll() {
1057        let e = create::<f64, f64>(
1058            ReferenceCellType::Quadrilateral,
1059            4,
1060            Continuity::Discontinuous,
1061            Variant::GLL,
1062        );
1063
1064        let mut data = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, 4));
1065        let mut points = rlst_dynamic_array!(f64, [2, 4]);
1066        *points.get_mut([0, 0]).unwrap() = 0.25;
1067        *points.get_mut([1, 0]).unwrap() = 0.25;
1068        *points.get_mut([0, 1]).unwrap() = 0.25;
1069        *points.get_mut([1, 1]).unwrap() = 0.75;
1070        *points.get_mut([0, 2]).unwrap() = 0.75;
1071        *points.get_mut([1, 2]).unwrap() = 0.25;
1072        *points.get_mut([0, 3]).unwrap() = 0.75;
1073        *points.get_mut([1, 3]).unwrap() = 0.75;
1074
1075        e.tabulate(&points, 0, &mut data);
1076
1077        for pt in 0..data.shape()[1] {
1078            for basis in 0..data.shape()[1] {
1079                assert!(data.get([0, pt, basis, 0]).unwrap().abs() > 1e-10);
1080            }
1081        }
1082    }
1083}
ndelement/ciarlet/lagrange.rs

ndelement/ciarlet/
lagrange.rs
