Struct RaviartThomasElementFamily 

pub struct RaviartThomasElementFamily<T: RlstScalar + Getrf + Getri = f64, TGeo: RlstScalar = f64> { /* private fields */ }

Raviart-Thomas element family

A family of Raviart-Thomas elements on multiple cell types with appropriate continuity across different cell types.

Implementations

impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> RaviartThomasElementFamily<T, TGeo>

pub fn new(degree: usize, continuity: Continuity) -> Self

Create new family with given degree and continuity.

Examples found in repository

ndfunctionspace/examples/test_mass_matrix.rs (line 90)

fn test_rt_mass_matrix() {
    let mesh = regular_sphere(0, ReferenceCellType::Triangle);

    let family = RaviartThomasElementFamily::<f64>::new(1, Continuity::Standard);
    let space = FunctionSpaceImpl::new(&mesh, &family);

    let mut mass_matrix = rlst_dynamic_array!(f64, [space.local_size(), space.local_size()]);

    let element = &space.elements()[0];

    let (p, w) = single_integral_quadrature(
        QuadratureRule::XiaoGimbutas,
        Domain::Triangle,
        2 * element.lagrange_superdegree(),
    )
    .unwrap();
    let npts = w.len();
    let mut pts = rlst_dynamic_array!(f64, [2, npts]);
    for i in 0..w.len() {
        for j in 0..2 {
            *pts.get_mut([j, i]).unwrap() = p[3 * i + j];
        }
    }
    let wts = w.iter().map(|i| *i / 2.0).collect::<Vec<_>>();

    let mut table = DynArray::<f64, 4>::from_shape(element.tabulate_array_shape(0, npts));
    element.tabulate(&pts, 0, &mut table);
    let mut pushed_table = rlst_dynamic_array!(
        f64,
        [table.shape()[0], table.shape()[1], table.shape()[2], 3]
    );

    let gmap = mesh.geometry_map(ReferenceCellType::Triangle, 1, &pts);
    let mut jacobians = rlst_dynamic_array!(f64, [mesh.geometry_dim(), mesh.topology_dim(), npts]);
    let mut jinv = rlst_dynamic_array!(f64, [mesh.topology_dim(), mesh.geometry_dim(), npts]);
    let mut jdets = vec![0.0; npts];

    for cell in mesh.entity_iter(ReferenceCellType::Triangle) {
        let dofs = space
            .entity_closure_dofs(ReferenceCellType::Triangle, cell.local_index())
            .unwrap();
        gmap.jacobians_inverses_dets(cell.local_index(), &mut jacobians, &mut jinv, &mut jdets);
        element.push_forward(&table, 0, &jacobians, &jdets, &jinv, &mut pushed_table);

        for (test_i, test_dof) in dofs.iter().enumerate() {
            for (trial_i, trial_dof) in dofs.iter().enumerate() {
                *mass_matrix.get_mut([*test_dof, *trial_dof]).unwrap() += wts
                    .iter()
                    .enumerate()
                    .map(|(i, w)| {
                        jdets[i]
                            * *w
                            * (0..3)
                                .map(|j| {
                                    *pushed_table.get([0, i, test_i, j]).unwrap()
                                        * *pushed_table.get([0, i, trial_i, j]).unwrap()
                                })
                                .sum::<f64>()
                    })
                    .sum::<f64>();
            }
        }
    }

    // Compare matrix entries to FEniCS
    for i in 0..12 {
        assert_relative_eq!(mass_matrix[[i, i]], 0.4811252243246884, epsilon = 1e-10);
    }
    for i in 0..12 {
        for j in 0..12 {
            if i != j && mass_matrix[[i, j]].abs() > 0.001 {
                assert_relative_eq!(
                    mass_matrix[[i, j]].abs(),
                    0.0481125224324689,
                    epsilon = 1e-10
                );
            }
        }
    }
}

Trait Implementations

impl<T: RlstScalar + Getrf + Getri, TGeo: RlstScalar> ElementFamily for RaviartThomasElementFamily<T, TGeo>

type T = T

The scalar type

type CellType = ReferenceCellType

Cell type

type FiniteElement = CiarletElement<T, ContravariantPiolaMap, TGeo>

The finite element type

fn element( &self, cell_type: ReferenceCellType, ) -> CiarletElement<T, ContravariantPiolaMap, TGeo>

Create an element for the given cell type.