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