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 rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri};
10use rlst::{RlstScalar, rlst_dynamic_array};
11use std::marker::PhantomData;
12
13/// Create a Lagrange element.
14pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
15    cell_type: ReferenceCellType,
16    degree: usize,
17    continuity: Continuity,
18) -> CiarletElement<T, IdentityMap, TGeo> {
19    let dim = polynomial_count(cell_type, degree);
20    let tdim = reference_cell::dim(cell_type);
21    let mut wcoeffs = rlst_dynamic_array!(T, [dim, 1, dim]);
22    for i in 0..dim {
23        *wcoeffs.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap();
24    }
25
26    let mut x = [vec![], vec![], vec![], vec![]];
27    let mut m = [vec![], vec![], vec![], vec![]];
28    let entity_counts = reference_cell::entity_counts(cell_type);
29    let vertices = reference_cell::vertices::<TGeo>(cell_type);
30    if degree == 0 {
31        if continuity == Continuity::Standard {
32            panic!("Cannot create continuous degree 0 Lagrange element");
33        }
34        for (d, counts) in entity_counts.iter().enumerate() {
35            for _e in 0..*counts {
36                x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
37                m[d].push(rlst_dynamic_array!(T, [0, 1, 0]));
38            }
39        }
40        let mut midp = rlst_dynamic_array!(TGeo, [tdim, 1]);
41        let nvertices = entity_counts[0];
42        for i in 0..tdim {
43            for vertex in &vertices {
44                *midp.get_mut([i, 0]).unwrap() += num::cast::<_, TGeo>(vertex[i]).unwrap();
45            }
46            *midp.get_mut([i, 0]).unwrap() /= num::cast::<_, TGeo>(nvertices).unwrap();
47        }
48        x[tdim].push(midp);
49        let mut mentry = rlst_dynamic_array!(T, [1, 1, 1]);
50        *mentry.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
51        m[tdim].push(mentry);
52    } else {
53        let edges = reference_cell::edges(cell_type);
54        let faces = reference_cell::faces(cell_type);
55        let volumes = reference_cell::volumes(cell_type);
56        for vertex in &vertices {
57            let mut pts = rlst_dynamic_array!(TGeo, [tdim, 1]);
58            for (i, v) in vertex.iter().enumerate() {
59                *pts.get_mut([i, 0]).unwrap() = num::cast::<_, TGeo>(*v).unwrap();
60            }
61            x[0].push(pts);
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[0].push(mentry);
65        }
66        for e in &edges {
67            let mut pts = rlst_dynamic_array!(TGeo, [tdim, degree - 1]);
68            let [vn0, vn1] = e[..] else {
69                panic!();
70            };
71            let v0 = &vertices[vn0];
72            let v1 = &vertices[vn1];
73            let mut ident = rlst_dynamic_array!(T, [degree - 1, 1, degree - 1]);
74
75            for i in 1..degree {
76                *ident.get_mut([i - 1, 0, i - 1]).unwrap() = T::from(1.0).unwrap();
77                for j in 0..tdim {
78                    *pts.get_mut([j, i - 1]).unwrap() = num::cast::<_, TGeo>(v0[j]).unwrap()
79                        + num::cast::<_, TGeo>(i).unwrap() / num::cast::<_, TGeo>(degree).unwrap()
80                            * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap();
81                }
82            }
83            x[1].push(pts);
84            m[1].push(ident);
85        }
86        for (e, face_type) in reference_cell::entity_types(cell_type)[2]
87            .iter()
88            .enumerate()
89        {
90            let npts = match face_type {
91                ReferenceCellType::Triangle => {
92                    if degree > 2 {
93                        (degree - 1) * (degree - 2) / 2
94                    } else {
95                        0
96                    }
97                }
98                ReferenceCellType::Quadrilateral => (degree - 1).pow(2),
99                _ => {
100                    panic!("Unsupported face type");
101                }
102            };
103            let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]);
104
105            let [vn0, vn1, vn2] = faces[e][..3] else {
106                panic!();
107            };
108            let v0 = &vertices[vn0];
109            let v1 = &vertices[vn1];
110            let v2 = &vertices[vn2];
111
112            match face_type {
113                ReferenceCellType::Triangle => {
114                    let mut n = 0;
115                    for i0 in 1..degree {
116                        for i1 in 1..degree - i0 {
117                            for j in 0..tdim {
118                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
119                                    .unwrap()
120                                    + num::cast::<_, TGeo>(i0).unwrap()
121                                        / num::cast::<_, TGeo>(degree).unwrap()
122                                        * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
123                                    + num::cast::<_, TGeo>(i1).unwrap()
124                                        / num::cast::<_, TGeo>(degree).unwrap()
125                                        * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap();
126                            }
127                            n += 1;
128                        }
129                    }
130                }
131                ReferenceCellType::Quadrilateral => {
132                    let mut n = 0;
133                    for i0 in 1..degree {
134                        for i1 in 1..degree {
135                            for j in 0..tdim {
136                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
137                                    .unwrap()
138                                    + num::cast::<_, TGeo>(i0).unwrap()
139                                        / num::cast::<_, TGeo>(degree).unwrap()
140                                        * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
141                                    + num::cast::<_, TGeo>(i1).unwrap()
142                                        / num::cast::<_, TGeo>(degree).unwrap()
143                                        * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap();
144                            }
145                            n += 1;
146                        }
147                    }
148                }
149                _ => {
150                    panic!("Unsupported face type.");
151                }
152            };
153
154            let mut ident = rlst_dynamic_array!(T, [npts, 1, npts]);
155            for i in 0..npts {
156                *ident.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap();
157            }
158            x[2].push(pts);
159            m[2].push(ident);
160        }
161        for (e, volume_type) in reference_cell::entity_types(cell_type)[3]
162            .iter()
163            .enumerate()
164        {
165            let npts = match volume_type {
166                ReferenceCellType::Tetrahedron => {
167                    if degree > 2 {
168                        (degree - 1) * (degree - 2) * (degree - 3) / 6
169                    } else {
170                        0
171                    }
172                }
173                ReferenceCellType::Hexahedron => (degree - 1).pow(3),
174                _ => {
175                    panic!("Unsupported face type");
176                }
177            };
178            let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]);
179
180            match volume_type {
181                ReferenceCellType::Tetrahedron => {
182                    let [vn0, vn1, vn2, vn3] = volumes[e][..4] else {
183                        panic!();
184                    };
185                    let v0 = &vertices[vn0];
186                    let v1 = &vertices[vn1];
187                    let v2 = &vertices[vn2];
188                    let v3 = &vertices[vn3];
189
190                    let mut n = 0;
191                    for i0 in 1..degree {
192                        for i1 in 1..degree - i0 {
193                            for i2 in 1..degree - i0 - i1 {
194                                for j in 0..tdim {
195                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
196                                        .unwrap()
197                                        + num::cast::<_, TGeo>(i0).unwrap()
198                                            / num::cast::<_, TGeo>(degree).unwrap()
199                                            * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
200                                        + num::cast::<_, TGeo>(i1).unwrap()
201                                            / num::cast::<_, TGeo>(degree).unwrap()
202                                            * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap()
203                                        + num::cast::<_, TGeo>(i2).unwrap()
204                                            / num::cast::<_, TGeo>(degree).unwrap()
205                                            * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap();
206                                }
207                                n += 1;
208                            }
209                        }
210                    }
211                }
212                ReferenceCellType::Hexahedron => {
213                    let [vn0, vn1, vn2, _, vn3] = volumes[e][..5] else {
214                        panic!();
215                    };
216                    let v0 = &vertices[vn0];
217                    let v1 = &vertices[vn1];
218                    let v2 = &vertices[vn2];
219                    let v3 = &vertices[vn3];
220
221                    let mut n = 0;
222                    for i0 in 1..degree {
223                        for i1 in 1..degree {
224                            for i2 in 1..degree {
225                                for j in 0..tdim {
226                                    *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j])
227                                        .unwrap()
228                                        + num::cast::<_, TGeo>(i0).unwrap()
229                                            / num::cast::<_, TGeo>(degree).unwrap()
230                                            * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap()
231                                        + num::cast::<_, TGeo>(i1).unwrap()
232                                            / num::cast::<_, TGeo>(degree).unwrap()
233                                            * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap()
234                                        + num::cast::<_, TGeo>(i2).unwrap()
235                                            / num::cast::<_, TGeo>(degree).unwrap()
236                                            * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap();
237                                }
238                                n += 1;
239                            }
240                        }
241                    }
242                }
243                _ => {
244                    panic!("Unsupported face type.");
245                }
246            };
247
248            let mut ident = rlst_dynamic_array!(T, [npts, 1, npts]);
249            for i in 0..npts {
250                *ident.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap();
251            }
252            x[3].push(pts);
253            m[3].push(ident);
254        }
255    }
256    CiarletElement::<T, IdentityMap, TGeo>::create(
257        "Lagrange".to_string(),
258        cell_type,
259        degree,
260        vec![],
261        wcoeffs,
262        x,
263        m,
264        continuity,
265        degree,
266        IdentityMap {},
267    )
268}
269
270/// Lagrange element family.
271///
272/// A family of Lagrange elements on multiple cell types with appropriate
273/// continuity across different cell types.
274pub struct LagrangeElementFamily<T: RlstScalar + Getrf + Getri = f64, TGeo: RlstScalar = f64> {
275    degree: usize,
276    continuity: Continuity,
277    _t: PhantomData<T>,
278    _tgeo: PhantomData<TGeo>,
279}
280
281impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> LagrangeElementFamily<T, TGeo> {
282    /// Create new family with given `degree` and `continuity`.
283    pub fn new(degree: usize, continuity: Continuity) -> Self {
284        Self {
285            degree,
286            continuity,
287            _t: PhantomData,
288            _tgeo: PhantomData,
289        }
290    }
291}
292
293impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> ElementFamily
294    for LagrangeElementFamily<T, TGeo>
295{
296    type T = T;
297    type FiniteElement = CiarletElement<T, IdentityMap, TGeo>;
298    type CellType = ReferenceCellType;
299    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, IdentityMap, TGeo> {
300        create::<T, TGeo>(cell_type, self.degree, self.continuity)
301    }
302}
ndelement/ciarlet/lagrange.rs

ndelement/ciarlet/
lagrange.rs
