ndelement/ciarlet/
nedelec.rs

1//! Nedelec elements.
2
3use super::CiarletElement;
4use crate::map::CovariantPiolaMap;
5use crate::math::orthogonalise_3;
6use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials};
7use crate::quadrature::gauss_jacobi_rule;
8use crate::reference_cell;
9use crate::traits::ElementFamily;
10use crate::types::{Continuity, ReferenceCellType};
11use itertools::izip;
12use rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri};
13use rlst::{DynArray, RlstScalar, SliceArray, rlst_dynamic_array};
14use std::marker::PhantomData;
15
16fn create_simplex<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
17    cell_type: ReferenceCellType,
18    degree: usize,
19    continuity: Continuity,
20) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
21    if cell_type != ReferenceCellType::Triangle && cell_type != ReferenceCellType::Tetrahedron {
22        panic!("Invalid cell: {cell_type:?}");
23    }
24
25    if degree < 1 {
26        panic!("Degree must be at least 1");
27    }
28
29    let tdim = reference_cell::dim(cell_type);
30    let cell_q = gauss_jacobi_rule(cell_type, 2 * degree).unwrap();
31    let pts_t = cell_q
32        .points
33        .iter()
34        .map(|i| TGeo::from(*i).unwrap())
35        .collect::<Vec<_>>();
36    let pts = SliceArray::<TGeo, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
37
38    let mut phi = DynArray::<T, 3>::from_shape(legendre_shape(cell_type, &pts, degree, 0));
39    tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
40
41    let pdim = phi.shape()[1];
42
43    let pdim_facet_minus1 = polynomial_count(
44        if tdim == 3 {
45            ReferenceCellType::Triangle
46        } else {
47            ReferenceCellType::Interval
48        },
49        degree - 1,
50    );
51    let pdim_edge_minus1 = polynomial_count(ReferenceCellType::Interval, degree - 1);
52    let pdim_face_minus2 = if degree < 2 {
53        0
54    } else {
55        polynomial_count(ReferenceCellType::Triangle, degree - 2)
56    };
57    let pdim_minus1 = polynomial_count(cell_type, degree - 1);
58    let pdim_minus2 = if degree < 2 {
59        0
60    } else {
61        polynomial_count(cell_type, degree - 2)
62    };
63    let pdim_minus3 = if degree < 3 {
64        0
65    } else {
66        polynomial_count(cell_type, degree - 3)
67    };
68
69    let wdim = pdim_minus1 * tdim + pdim_facet_minus1 * (tdim - 1) + pdim_edge_minus1 * (tdim - 2);
70    let mut wcoeffs = rlst_dynamic_array!(T, [wdim, tdim, pdim]);
71
72    // vector polynomials of degree <= n-1
73    for i in 0..tdim {
74        for j in 0..pdim_minus1 {
75            *wcoeffs.get_mut([i * pdim_minus1 + j, i, j]).unwrap() = T::from(1.0).unwrap();
76        }
77    }
78
79    if cell_type == ReferenceCellType::Triangle {
80        // Additional Nedelec polynomials
81        for i in 0..pdim_facet_minus1 {
82            for j in pdim_minus1..pdim {
83                let w0 = cell_q
84                    .weights
85                    .iter()
86                    .enumerate()
87                    .map(|(w_i, wt)| {
88                        T::from(*wt).unwrap()
89                            * phi[[0, pdim_minus2 + i, w_i]]
90                            * T::from(pts[[0, w_i]]).unwrap()
91                            * phi[[0, j, w_i]]
92                    })
93                    .sum::<T>();
94                let w1 = cell_q
95                    .weights
96                    .iter()
97                    .enumerate()
98                    .map(|(w_i, wt)| {
99                        T::from(*wt).unwrap()
100                            * phi[[0, pdim_minus2 + i, w_i]]
101                            * T::from(pts[[1, w_i]]).unwrap()
102                            * phi[[0, j, w_i]]
103                    })
104                    .sum::<T>();
105                wcoeffs[[2 * pdim_minus1 + i, 0, j]] = w1;
106                wcoeffs[[2 * pdim_minus1 + i, 1, j]] = -w0;
107            }
108        }
109    } else {
110        // Additional Nedelec polynomials
111        for i in 0..pdim_facet_minus1 {
112            for j in pdim_minus1..pdim {
113                let w0 = cell_q
114                    .weights
115                    .iter()
116                    .enumerate()
117                    .map(|(w_i, wt)| {
118                        T::from(*wt).unwrap()
119                            * phi[[0, pdim_minus2 + i, w_i]]
120                            * T::from(pts[[0, w_i]]).unwrap()
121                            * phi[[0, j, w_i]]
122                    })
123                    .sum::<T>();
124                let w1 = cell_q
125                    .weights
126                    .iter()
127                    .enumerate()
128                    .map(|(w_i, wt)| {
129                        T::from(*wt).unwrap()
130                            * phi[[0, pdim_minus2 + i, w_i]]
131                            * T::from(pts[[1, w_i]]).unwrap()
132                            * phi[[0, j, w_i]]
133                    })
134                    .sum::<T>();
135                let w2 = cell_q
136                    .weights
137                    .iter()
138                    .enumerate()
139                    .map(|(w_i, wt)| {
140                        T::from(*wt).unwrap()
141                            * phi[[0, pdim_minus2 + i, w_i]]
142                            * T::from(pts[[2, w_i]]).unwrap()
143                            * phi[[0, j, w_i]]
144                    })
145                    .sum::<T>();
146
147                if i >= pdim_face_minus2 {
148                    wcoeffs[[tdim * pdim_minus1 + i - pdim_face_minus2, 2, j]] = w1;
149                    wcoeffs[[tdim * pdim_minus1 + i - pdim_face_minus2, 1, j]] = -w2;
150                }
151                wcoeffs[[
152                    tdim * pdim_minus1 + i + pdim_facet_minus1 - pdim_face_minus2,
153                    0,
154                    j,
155                ]] = w2;
156                wcoeffs[[
157                    tdim * pdim_minus1 + i + pdim_facet_minus1 * 2 - pdim_face_minus2,
158                    0,
159                    j,
160                ]] = -w1;
161                wcoeffs[[
162                    tdim * pdim_minus1 + i + pdim_facet_minus1 - pdim_face_minus2,
163                    2,
164                    j,
165                ]] = -w0;
166                wcoeffs[[
167                    tdim * pdim_minus1 + i + pdim_facet_minus1 * 2 - pdim_face_minus2,
168                    1,
169                    j,
170                ]] = w0;
171            }
172        }
173    };
174
175    orthogonalise_3(&mut wcoeffs, pdim_minus1 * tdim);
176
177    let mut x = [vec![], vec![], vec![], vec![]];
178    let mut m = [vec![], vec![], vec![], vec![]];
179
180    let entity_counts = reference_cell::entity_counts(cell_type);
181    let vertices = reference_cell::vertices::<TGeo>(cell_type);
182
183    for _ in 0..entity_counts[0] {
184        x[0].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
185        m[0].push(rlst_dynamic_array!(T, [0, tdim, 0]));
186    }
187
188    // DOFs on edges
189    let edge_q = gauss_jacobi_rule(ReferenceCellType::Interval, 2 * degree - 1).unwrap();
190    let edge_pts_t = edge_q
191        .points
192        .iter()
193        .map(|i| TGeo::from(*i).unwrap())
194        .collect::<Vec<_>>();
195    let edge_pts = SliceArray::<TGeo, 2>::from_shape(&edge_pts_t, [1, edge_q.npoints]);
196
197    let mut edge_phi = DynArray::<T, 3>::from_shape(legendre_shape(
198        ReferenceCellType::Interval,
199        &edge_pts,
200        degree - 1,
201        0,
202    ));
203    tabulate_legendre_polynomials(
204        ReferenceCellType::Interval,
205        &edge_pts,
206        degree - 1,
207        0,
208        &mut edge_phi,
209    );
210
211    for edge in reference_cell::edges(cell_type) {
212        let mut pts = rlst_dynamic_array!(TGeo, [tdim, edge_q.npoints]);
213        let mut mat = rlst_dynamic_array!(T, [pdim_edge_minus1, tdim, edge_q.npoints]);
214
215        for (w_i, (pt, wt)) in izip!(&edge_pts_t, &edge_q.weights).enumerate() {
216            for i in 0..tdim {
217                pts[[i, w_i]] =
218                    vertices[edge[0]][i] + (vertices[edge[1]][i] - vertices[edge[0]][i]) * *pt;
219
220                for j in 0..pdim_edge_minus1 {
221                    mat[[j, i, w_i]] = T::from(*wt).unwrap()
222                        * edge_phi[[0, j, w_i]]
223                        * T::from(vertices[edge[1]][i] - vertices[edge[0]][i]).unwrap();
224                }
225            }
226        }
227
228        x[1].push(pts);
229        m[1].push(mat);
230    }
231
232    // DOFs on faces
233    if degree == 1 {
234        for _ in 0..entity_counts[2] {
235            x[2].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
236            m[2].push(rlst_dynamic_array!(T, [0, tdim, 0]))
237        }
238    } else {
239        let face_q = gauss_jacobi_rule(ReferenceCellType::Triangle, 2 * degree - 2).unwrap();
240        let face_pts_t = face_q
241            .points
242            .iter()
243            .map(|i| TGeo::from(*i).unwrap())
244            .collect::<Vec<_>>();
245        let face_pts = SliceArray::<TGeo, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
246
247        let mut face_phi = DynArray::<T, 3>::from_shape(legendre_shape(
248            ReferenceCellType::Triangle,
249            &face_pts,
250            degree - 2,
251            0,
252        ));
253        tabulate_legendre_polynomials(
254            ReferenceCellType::Triangle,
255            &face_pts,
256            degree - 2,
257            0,
258            &mut face_phi,
259        );
260
261        for face in reference_cell::faces(cell_type) {
262            let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]);
263            let mut mat = rlst_dynamic_array!(T, [2 * pdim_face_minus2, tdim, face_q.npoints]);
264
265            for (w_i, wt) in face_q.weights.iter().enumerate() {
266                for i in 0..tdim {
267                    pts[[i, w_i]] = vertices[face[0]][i]
268                        + (vertices[face[1]][i] - vertices[face[0]][i]) * face_pts[[0, w_i]]
269                        + (vertices[face[2]][i] - vertices[face[0]][i]) * face_pts[[1, w_i]];
270
271                    for tangent in 0..2 {
272                        for j in 0..pdim_face_minus2 {
273                            mat[[tangent * pdim_face_minus2 + j, i, w_i]] = T::from(*wt).unwrap()
274                                * face_phi[[0, j, w_i]]
275                                * T::from(vertices[face[tangent + 1]][i] - vertices[face[0]][i])
276                                    .unwrap();
277                        }
278                    }
279                }
280            }
281            x[2].push(pts);
282            m[2].push(mat);
283        }
284    }
285    if tdim == 3 {
286        if degree <= 2 {
287            x[3].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
288            m[3].push(rlst_dynamic_array!(T, [0, tdim, 0]))
289        } else {
290            let internal_q = gauss_jacobi_rule(cell_type, 2 * degree - 3).unwrap();
291            let mut pts = rlst_dynamic_array!(TGeo, [tdim, internal_q.npoints]);
292            for p in 0..internal_q.npoints {
293                for d in 0..3 {
294                    pts[[d, p]] = TGeo::from(internal_q.points[3 * p + d]).unwrap()
295                }
296            }
297
298            let mut internal_phi =
299                DynArray::<T, 3>::from_shape(legendre_shape(cell_type, &pts, degree - 3, 0));
300            tabulate_legendre_polynomials(cell_type, &pts, degree - 3, 0, &mut internal_phi);
301
302            let mut mat = rlst_dynamic_array!(T, [tdim * pdim_minus3, tdim, internal_q.npoints]);
303            for (w_i, wt) in internal_q.weights.iter().enumerate() {
304                for i in 0..tdim {
305                    for j in 0..pdim_minus3 {
306                        mat[[j + pdim_minus3 * i, i, w_i]] =
307                            T::from(*wt).unwrap() * internal_phi[[0, j, w_i]];
308                    }
309                }
310            }
311
312            x[tdim].push(pts);
313            m[tdim].push(mat);
314        }
315    }
316
317    CiarletElement::create(
318        "Nedelec (first kind)".to_string(),
319        cell_type,
320        degree,
321        vec![tdim],
322        wcoeffs,
323        x,
324        m,
325        continuity,
326        degree,
327        CovariantPiolaMap {},
328    )
329}
330
331fn create_tp<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
332    cell_type: ReferenceCellType,
333    degree: usize,
334    continuity: Continuity,
335) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
336    if cell_type != ReferenceCellType::Quadrilateral && cell_type != ReferenceCellType::Hexahedron {
337        panic!("Invalid cell: {cell_type:?}");
338    }
339
340    if degree < 1 {
341        panic!("Degree must be at least 1");
342    }
343
344    let tdim = reference_cell::dim(cell_type);
345    let cell_q = gauss_jacobi_rule(cell_type, 2 * degree).unwrap();
346    let pts_t = cell_q
347        .points
348        .iter()
349        .map(|i| TGeo::from(*i).unwrap())
350        .collect::<Vec<_>>();
351    let pts = SliceArray::<TGeo, 2>::from_shape(&pts_t, [tdim, cell_q.npoints]);
352
353    let mut phi = DynArray::<T, 3>::from_shape(legendre_shape(cell_type, &pts, degree, 0));
354    tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
355
356    let pdim = phi.shape()[1];
357
358    let pdim_edge = polynomial_count(ReferenceCellType::Interval, degree);
359    let pdim_edge_minus1 = polynomial_count(ReferenceCellType::Interval, degree - 1);
360    let pdim_edge_minus2 = if degree < 2 {
361        0
362    } else {
363        polynomial_count(ReferenceCellType::Interval, degree - 2)
364    };
365
366    let entity_counts = reference_cell::entity_counts(cell_type);
367
368    let wdim = entity_counts[1] * pdim_edge_minus1
369        + entity_counts[2] * 2 * pdim_edge_minus1 * pdim_edge_minus2
370        + entity_counts[3] * 3 * pdim_edge_minus1 * pdim_edge_minus2 * pdim_edge_minus2;
371    let mut wcoeffs = rlst_dynamic_array!(T, [wdim, tdim, pdim]);
372
373    // vector polynomials of degree <= n-1
374    if tdim == 2 {
375        for i in 0..pdim_edge_minus1 {
376            for j in 0..pdim_edge {
377                *wcoeffs
378                    .get_mut([i * pdim_edge + j, 0, i * pdim_edge + j])
379                    .unwrap() = T::from(1.0).unwrap();
380                *wcoeffs
381                    .get_mut([
382                        pdim_edge_minus1 * pdim_edge + j * pdim_edge_minus1 + i,
383                        1,
384                        j * pdim_edge + i,
385                    ])
386                    .unwrap() = T::from(1.0).unwrap();
387            }
388        }
389    } else {
390        for i in 0..pdim_edge_minus1 {
391            for j in 0..pdim_edge {
392                for k in 0..pdim_edge {
393                    *wcoeffs
394                        .get_mut([
395                            i * pdim_edge.pow(2) + j * pdim_edge + k,
396                            0,
397                            i * pdim_edge.pow(2) + j * pdim_edge + k,
398                        ])
399                        .unwrap() = T::from(1.0).unwrap();
400                    *wcoeffs
401                        .get_mut([
402                            pdim_edge.pow(2) * pdim_edge_minus1
403                                + k * pdim_edge * pdim_edge_minus1
404                                + i * pdim_edge
405                                + j,
406                            1,
407                            k * pdim_edge.pow(2) + i * pdim_edge + j,
408                        ])
409                        .unwrap() = T::from(1.0).unwrap();
410                    *wcoeffs
411                        .get_mut([
412                            pdim_edge.pow(2) * pdim_edge_minus1 * 2
413                                + j * pdim_edge * pdim_edge_minus1
414                                + k * pdim_edge_minus1
415                                + i,
416                            2,
417                            j * pdim_edge.pow(2) + k * pdim_edge + i,
418                        ])
419                        .unwrap() = T::from(1.0).unwrap();
420                }
421            }
422        }
423    }
424
425    let mut x = [vec![], vec![], vec![], vec![]];
426    let mut m = [vec![], vec![], vec![], vec![]];
427
428    let vertices = reference_cell::vertices::<TGeo>(cell_type);
429
430    for _ in 0..entity_counts[0] {
431        x[0].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
432        m[0].push(rlst_dynamic_array!(T, [0, tdim, 0]));
433    }
434
435    // DOFs on edges
436    let edge_q = gauss_jacobi_rule(ReferenceCellType::Interval, 2 * degree - 1).unwrap();
437    let edge_pts_t = edge_q
438        .points
439        .iter()
440        .map(|i| TGeo::from(*i).unwrap())
441        .collect::<Vec<_>>();
442    let edge_pts = SliceArray::<TGeo, 2>::from_shape(&edge_pts_t, [1, edge_q.npoints]);
443
444    let mut edge_phi = DynArray::<T, 3>::from_shape(legendre_shape(
445        ReferenceCellType::Interval,
446        &edge_pts,
447        degree - 1,
448        0,
449    ));
450    tabulate_legendre_polynomials(
451        ReferenceCellType::Interval,
452        &edge_pts,
453        degree - 1,
454        0,
455        &mut edge_phi,
456    );
457
458    for edge in reference_cell::edges(cell_type) {
459        let mut pts = rlst_dynamic_array!(TGeo, [tdim, edge_q.npoints]);
460        let mut mat = rlst_dynamic_array!(T, [pdim_edge_minus1, tdim, edge_q.npoints]);
461
462        for (w_i, (pt, wt)) in izip!(&edge_pts_t, &edge_q.weights).enumerate() {
463            for i in 0..tdim {
464                pts[[i, w_i]] =
465                    vertices[edge[0]][i] + (vertices[edge[1]][i] - vertices[edge[0]][i]) * *pt;
466
467                for j in 0..pdim_edge_minus1 {
468                    mat[[j, i, w_i]] = T::from(*wt).unwrap()
469                        * edge_phi[[0, j, w_i]]
470                        * T::from(vertices[edge[1]][i] - vertices[edge[0]][i]).unwrap();
471                }
472            }
473        }
474
475        x[1].push(pts);
476        m[1].push(mat);
477    }
478
479    // DOFs on faces
480    if degree == 1 {
481        for _ in 0..entity_counts[2] {
482            x[2].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
483            m[2].push(rlst_dynamic_array!(T, [0, tdim, 0]))
484        }
485    } else {
486        let face_q = gauss_jacobi_rule(ReferenceCellType::Quadrilateral, 2 * degree - 1).unwrap();
487        let face_pts_t = face_q
488            .points
489            .iter()
490            .map(|i| TGeo::from(*i).unwrap())
491            .collect::<Vec<_>>();
492        let face_pts = SliceArray::<TGeo, 2>::from_shape(&face_pts_t, [2, face_q.npoints]);
493
494        let mut face_phi = DynArray::<T, 3>::from_shape(legendre_shape(
495            ReferenceCellType::Quadrilateral,
496            &face_pts,
497            degree - 1,
498            0,
499        ));
500        tabulate_legendre_polynomials(
501            ReferenceCellType::Quadrilateral,
502            &face_pts,
503            degree - 1,
504            0,
505            &mut face_phi,
506        );
507
508        for face in reference_cell::faces(cell_type) {
509            let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]);
510            let mdim = 2 * pdim_edge_minus2 * pdim_edge_minus1;
511            let mut mat = rlst_dynamic_array!(T, [mdim, tdim, face_q.npoints]);
512
513            for (w_i, wt) in face_q.weights.iter().enumerate() {
514                for i in 0..tdim {
515                    pts[[i, w_i]] = vertices[face[0]][i]
516                        + (vertices[face[1]][i] - vertices[face[0]][i]) * face_pts[[0, w_i]]
517                        + (vertices[face[2]][i] - vertices[face[0]][i]) * face_pts[[1, w_i]];
518                }
519                for i in 0..pdim_edge_minus2 {
520                    for j in 0..pdim_edge_minus1 {
521                        let index = 2 * (i * pdim_edge_minus1 + j);
522                        let entry0 =
523                            T::from(*wt).unwrap() * face_phi[[0, j * pdim_edge_minus1 + i, w_i]];
524                        let entry1 =
525                            T::from(*wt).unwrap() * face_phi[[0, i * pdim_edge_minus1 + j, w_i]];
526                        for d in 0..tdim {
527                            mat[[index, d, w_i]] = entry0
528                                * T::from(vertices[face[1]][d] - vertices[face[0]][d]).unwrap();
529                            mat[[index + 1, d, w_i]] = entry1
530                                * T::from(vertices[face[2]][d] - vertices[face[0]][d]).unwrap();
531                        }
532                    }
533                }
534            }
535            x[2].push(pts);
536            m[2].push(mat);
537        }
538    }
539    // DOFs on volume
540    if tdim == 3 {
541        if degree == 1 {
542            x[3].push(rlst_dynamic_array!(TGeo, [tdim, 0]));
543            m[3].push(rlst_dynamic_array!(T, [0, tdim, 0]))
544        } else {
545            let interior_q =
546                gauss_jacobi_rule(ReferenceCellType::Hexahedron, 2 * degree - 1).unwrap();
547            let interior_pts_t = interior_q
548                .points
549                .iter()
550                .map(|i| TGeo::from(*i).unwrap())
551                .collect::<Vec<_>>();
552            let interior_pts =
553                SliceArray::<TGeo, 2>::from_shape(&interior_pts_t, [3, interior_q.npoints]);
554
555            let mut interior_phi = DynArray::<T, 3>::from_shape(legendre_shape(
556                ReferenceCellType::Hexahedron,
557                &interior_pts,
558                degree - 1,
559                0,
560            ));
561            tabulate_legendre_polynomials(
562                ReferenceCellType::Hexahedron,
563                &interior_pts,
564                degree - 1,
565                0,
566                &mut interior_phi,
567            );
568
569            let mut pts = rlst_dynamic_array!(TGeo, [tdim, interior_q.npoints]);
570            let mdim = 3 * pdim_edge_minus2.pow(2) * pdim_edge_minus1;
571            let mut mat = rlst_dynamic_array!(T, [mdim, tdim, interior_q.npoints]);
572
573            for (w_i, wt) in interior_q.weights.iter().enumerate() {
574                for i in 0..tdim {
575                    pts[[i, w_i]] = interior_pts[[i, w_i]];
576                }
577                for i in 0..pdim_edge_minus2 {
578                    for j in 0..pdim_edge_minus2 {
579                        for k in 0..pdim_edge_minus1 {
580                            let index = 3
581                                * (i * pdim_edge_minus1 * pdim_edge_minus2
582                                    + j * pdim_edge_minus1
583                                    + k);
584                            mat[[index, 0, w_i]] = T::from(*wt).unwrap()
585                                * interior_phi[[
586                                    0,
587                                    k * pdim_edge_minus1.pow(2) + j * pdim_edge_minus1 + i,
588                                    w_i,
589                                ]];
590                            mat[[index + 1, 1, w_i]] = T::from(*wt).unwrap()
591                                * interior_phi[[
592                                    0,
593                                    i * pdim_edge_minus1.pow(2) + k * pdim_edge_minus1 + j,
594                                    w_i,
595                                ]];
596                            mat[[index + 2, 2, w_i]] = T::from(*wt).unwrap()
597                                * interior_phi[[
598                                    0,
599                                    j * pdim_edge_minus1.pow(2) + i * pdim_edge_minus1 + k,
600                                    w_i,
601                                ]];
602                        }
603                    }
604                }
605            }
606            x[3].push(pts);
607            m[3].push(mat);
608        }
609    }
610
611    CiarletElement::create(
612        "Nedelec (first kind)".to_string(),
613        cell_type,
614        degree,
615        vec![tdim],
616        wcoeffs,
617        x,
618        m,
619        continuity,
620        degree,
621        CovariantPiolaMap {},
622    )
623}
624
625/// Create a Nedelec (first kind) element
626pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
627    cell_type: ReferenceCellType,
628    degree: usize,
629    continuity: Continuity,
630) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
631    if cell_type == ReferenceCellType::Triangle || cell_type == ReferenceCellType::Tetrahedron {
632        create_simplex::<T, TGeo>(cell_type, degree, continuity)
633    } else if cell_type == ReferenceCellType::Quadrilateral
634        || cell_type == ReferenceCellType::Hexahedron
635    {
636        create_tp::<T, TGeo>(cell_type, degree, continuity)
637    } else {
638        panic!("Invalid cell: {cell_type:?}");
639    }
640}
641
642/// Nedelec (first kind) element family
643///
644/// A family of Nedelec elements on multiple cell types with appropriate
645/// continuity across different cell types.
646pub struct NedelecFirstKindElementFamily<
647    T: RlstScalar + Getrf + Getri = f64,
648    TGeo: RlstScalar = f64,
649> {
650    degree: usize,
651    continuity: Continuity,
652    _t: PhantomData<T>,
653    _tgeo: PhantomData<TGeo>,
654}
655
656impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> NedelecFirstKindElementFamily<T, TGeo> {
657    /// Create new family with given `degree` and `continuity`.
658    pub fn new(degree: usize, continuity: Continuity) -> Self {
659        Self {
660            degree,
661            continuity,
662            _t: PhantomData,
663            _tgeo: PhantomData,
664        }
665    }
666}
667
668impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> ElementFamily
669    for NedelecFirstKindElementFamily<T, TGeo>
670{
671    type T = T;
672    type CellType = ReferenceCellType;
673    type FiniteElement = CiarletElement<T, CovariantPiolaMap, TGeo>;
674    fn element(&self, cell_type: ReferenceCellType) -> CiarletElement<T, CovariantPiolaMap, TGeo> {
675        create::<T, TGeo>(cell_type, self.degree, self.continuity)
676    }
677}
ndelement/ciarlet/nedelec.rs

ndelement/ciarlet/
nedelec.rs
