Struct FunctionSpaceImpl 

pub struct FunctionSpaceImpl<'a, T: RlstScalar, TMesh: RlstScalar, E: Debug + PartialEq + Eq + Clone + Copy + Hash, G: Mesh<EntityDescriptor = E, T = TMesh>, F: MappedFiniteElement<CellType = E, T = T>> { /* private fields */ }

Function space.

Implementations

impl<'a, T: RlstScalar, TMesh: RlstScalar, G: Mesh<EntityDescriptor = ReferenceCellType, T = TMesh>, F: MappedFiniteElement<CellType = ReferenceCellType, T = T>> FunctionSpaceImpl<'a, T, TMesh, ReferenceCellType, G, F>

pub fn new<EF: ElementFamily<FiniteElement = F, CellType = ReferenceCellType>>( mesh: &'a G, family: &EF, ) -> Self

Create a new serial function space

Examples found in repository

ndfunctionspace/examples/test_parallel_space.rs (line 30)

fn test_parallel_function_space<C: Communicator>(comm: &C) {
    let mesh = unit_cube_distributed::<f64, _>(
        comm,
        GraphPartitioner::None,
        4,
        4,
        4,
        ReferenceCellType::Tetrahedron,
    );

    let family =
        LagrangeElementFamily::<f64>::new(2, Continuity::Standard, LagrangeVariant::Equispaced);
    let space = ParallelFunctionSpaceImpl::new(&mesh, &family);
    let serial_mesh = unit_cube::<f64>(4, 4, 4, ReferenceCellType::Tetrahedron);
    let serial_space = FunctionSpaceImpl::new(&serial_mesh, &family);

    assert_eq!(space.global_size(), serial_space.global_size());
}

ndfunctionspace/examples/test_mass_matrix.rs (line 23)

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

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

    let mut mass_matrix = rlst_dynamic_array!(f64, [space.process_size(), space.process_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 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);
        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
                            * *table.get([0, i, test_i, 0]).unwrap()
                            * *table.get([0, i, trial_i, 0]).unwrap()
                    })
                    .sum::<f64>();
            }
        }
    }

    // Compare matrix entries to values from Bempp
    for i in 0..6 {
        assert_relative_eq!(mass_matrix[[i, i]], 0.5773502691896255, epsilon = 1e-10);
    }
    for i in 0..6 {
        for j in 0..6 {
            if i != j && mass_matrix[[i, j]].abs() > 0.001 {
                assert_relative_eq!(mass_matrix[[i, j]], 0.1443375672974061, epsilon = 1e-10);
            }
        }
    }
}

/// Test values in Raviart-Thomas mass matrix
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.process_size(), space.process_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<'a, T: RlstScalar, TMesh: RlstScalar, E: Debug + PartialEq + Eq + Clone + Copy + Hash, G: Mesh<EntityDescriptor = E, T = TMesh>, F: MappedFiniteElement<CellType = E, T = T>> FunctionSpace for FunctionSpaceImpl<'a, T, TMesh, E, G, F>

type T = T

Scalar type

type TMesh = TMesh

Scalar type for geometry

type EntityDescriptor = E

Type used as identifier of different entity types

type Mesh = G

The type for the mesh this function space is defined on

type FiniteElement = F

The type for the finite element this mesh is defined by

fn mesh(&self) -> &G

The mesh that this function space is defined on

fn elements(&self) -> &[F]

The finite elements used in this function space

fn entities_by_element(&self, element_index: usize) -> Option<&[usize]>

A list of entity indices that use the element with the given index

fn entity_dofs(&self, entity_type: E, entity_number: usize) -> Option<&[usize]>

Get the process DOF numbers associated with the given entity

fn entity_closure_dofs( &self, entity_type: E, entity_number: usize, ) -> Option<&[usize]>

Get the process DOF numbers associated with the closure of the given entity

fn process_size(&self) -> usize

Get the number of process DOFs

fn process_owned_size(&self) -> usize

Get the number of owned process DOFs

fn global_size(&self) -> usize

Get the number of DOFs on all processes

fn global_dof_index(&self, process_dof_index: usize) -> usize

Get the global DOF index associated with a process DOF index

fn ownership(&self, _process_dof_index: usize) -> Ownership

Get ownership of a process DOF