pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
cell_type: ReferenceCellType,
degree: usize,
continuity: Continuity,
) -> CiarletElement<T, IdentityMap, TGeo>Expand description
Create a Lagrange element.
Examples found in repository?
ndelement/examples/lagrange_element.rs (line 11)
8fn main() {
9 // Create a P2 element on a triangle
10 let element =
11 lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
12
13 println!("This element has {} basis functions.", element.dim());
14
15 // Create an array to store the basis function values
16 let mut basis_values = DynArray::<f64, 4>::from_shape(element.tabulate_array_shape(0, 1));
17 // Create array containing the point [1/3, 1/3]
18 let mut points = rlst_dynamic_array!(f64, [2, 1]);
19 points[[0, 0]] = 1.0 / 3.0;
20 points[[1, 0]] = 1.0 / 3.0;
21 // Tabulate the element's basis functions at the point
22 element.tabulate(&points, 0, &mut basis_values);
23 println!(
24 "The values of the basis functions at the point (1/3, 1/3) are: {:?}",
25 basis_values.data()
26 );
27
28 // Set point to [1, 0]
29 points[[0, 0]] = 1.0;
30 points[[1, 0]] = 0.0;
31 // Tabulate the element's basis functions at the point
32 element.tabulate(&points, 0, &mut basis_values);
33 println!(
34 "The values of the basis functions at the point (1, 0) are: {:?}",
35 basis_values.data()
36 );
37}More examples
ndgrid/examples/test_push_forward.rs (line 25)
15fn test_lagrange_push_forward() {
16 let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
17 b.add_point(0, &[0.0, 0.0, 0.0]);
18 b.add_point(1, &[1.0, 0.0, 0.0]);
19 b.add_point(2, &[2.0, 0.0, 1.0]);
20 b.add_point(3, &[0.0, 1.0, 0.0]);
21 b.add_cell(0, &[0, 1, 3]);
22 b.add_cell(1, &[1, 2, 3]);
23 let grid = b.create_grid();
24
25 let e = lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 4, Continuity::Standard);
26
27 let npts = 5;
28
29 let mut cell0_points = rlst_dynamic_array!(f64, [2, npts]);
30 for i in 0..npts {
31 cell0_points[[1, i]] = i as f64 / (npts - 1) as f64;
32 cell0_points[[0, i]] = 1.0 - cell0_points[[1, i]];
33 }
34 let mut cell0_table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
35 e.tabulate(&cell0_points, 0, &mut cell0_table);
36
37 let mut points = rlst_dynamic_array!(f64, [2, npts]);
38 for i in 0..npts {
39 points[[0, i]] = 0.0;
40 points[[1, i]] = i as f64 / (npts - 1) as f64;
41 }
42
43 let mut table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
44 e.tabulate(&points, 0, &mut table);
45
46 let gmap = grid.geometry_map(ReferenceCellType::Triangle, 1, &points);
47
48 let mut jacobians = rlst_dynamic_array!(f64, [grid.geometry_dim(), grid.topology_dim(), npts]);
49 let mut jinv = rlst_dynamic_array!(f64, [grid.topology_dim(), grid.geometry_dim(), npts]);
50 let mut jdets = vec![0.0; npts];
51
52 gmap.jacobians_inverses_dets(1, &mut jacobians, &mut jinv, &mut jdets);
53
54 let mut cell1_table = DynArray::<f64, 4>::from_shape(table.shape());
55 e.push_forward(&table, 0, &jacobians, &jdets, &jinv, &mut cell1_table);
56
57 // Check that basis functions are continuous between cells
58 for (cell0_dof, cell1_dof) in izip!(
59 e.entity_closure_dofs(1, 1).unwrap(),
60 e.entity_closure_dofs(1, 0).unwrap()
61 ) {
62 for i in 0..npts {
63 assert_relative_eq!(
64 cell0_table[[0, i, *cell0_dof, 0]],
65 cell1_table[[0, i, *cell1_dof, 0]],
66 epsilon = 1e-10
67 );
68 }
69 }
70}