pub fn create<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar>(
cell_type: ReferenceCellType,
degree: usize,
continuity: Continuity,
) -> CiarletElement<T, ContravariantPiolaMap, TGeo>Expand description
Create a Raviart-Thomas element
Examples found in repository?
ndmesh/examples/test_push_forward.rs (line 89)
78fn test_rt_push_forward() {
79 let mut b = SingleElementMeshBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
80 b.add_point(0, &[0.0, 0.0, 0.0]);
81 b.add_point(1, &[1.0, 0.0, 0.0]);
82 b.add_point(2, &[2.0, 0.0, 1.0]);
83 b.add_point(3, &[0.0, 1.0, 0.0]);
84 b.add_cell(0, &[0, 1, 3]);
85 b.add_cell(1, &[1, 2, 3]);
86 let mesh = b.create_mesh();
87
88 let e =
89 raviart_thomas::create::<f64, f64>(ReferenceCellType::Triangle, 4, Continuity::Standard);
90
91 let npts = 5;
92
93 let mut cell0_points = rlst_dynamic_array!(f64, [2, npts]);
94 for i in 0..npts {
95 cell0_points[[1, i]] = i as f64 / (npts - 1) as f64;
96 cell0_points[[0, i]] = 1.0 - cell0_points[[1, i]];
97 }
98 let mut cell0_table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
99 e.tabulate(&cell0_points, 0, &mut cell0_table);
100
101 let mut points = rlst_dynamic_array!(f64, [2, npts]);
102 for i in 0..npts {
103 points[[0, i]] = 0.0;
104 points[[1, i]] = i as f64 / (npts - 1) as f64;
105 }
106
107 let mut table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
108 e.tabulate(&points, 0, &mut table);
109
110 let gmap = mesh.geometry_map(ReferenceCellType::Triangle, 1, &points);
111
112 let mut jacobians = rlst_dynamic_array!(f64, [mesh.geometry_dim(), mesh.topology_dim(), npts]);
113 let mut jinv = rlst_dynamic_array!(f64, [mesh.topology_dim(), mesh.geometry_dim(), npts]);
114 let mut jdets = vec![0.0; npts];
115
116 gmap.jacobians_inverses_dets(1, &mut jacobians, &mut jinv, &mut jdets);
117
118 let mut cell1_table =
119 DynArray::<f64, 4>::from_shape([table.shape()[0], table.shape()[1], table.shape()[2], 3]);
120 e.push_forward(&table, 0, &jacobians, &jdets, &jinv, &mut cell1_table);
121
122 // Check that basis functions dotted with normal to edge are continuous between cells
123 for (cell0_dof, cell1_dof) in izip!(
124 e.entity_closure_dofs(1, 0).unwrap(),
125 e.entity_closure_dofs(1, 1).unwrap()
126 ) {
127 for i in 0..npts {
128 assert_relative_eq!(
129 (cell0_table[[0, i, *cell0_dof, 0]] + cell0_table[[0, i, *cell0_dof, 1]])
130 / f64::sqrt(2.0),
131 (cell1_table[[0, i, *cell1_dof, 0]]
132 + cell1_table[[0, i, *cell1_dof, 1]]
133 + 2.0 * cell1_table[[0, i, *cell1_dof, 2]])
134 / f64::sqrt(6.0),
135 epsilon = 1e-10
136 );
137 }
138 }
139}