Struct SingleElementGridBuilder 

pub struct SingleElementGridBuilder<T: Scalar> { /* private fields */ }

Grid builder for a single element grid

The following gives an example of creating a new grid consisting of a single triangle.

use ndgrid::traits::Builder;
use ndgrid::SingleElementGridBuilder;
use ndelement::types::ReferenceCellType;

// The geometric dimension of our space is 3.
let gdim = 3;

// We are building a two dimensional surface triangle grid within a three dimensional space.
// Our grid will have three points and one `Triangle` cell of order 1.
let mut builder = SingleElementGridBuilder::new_with_capacity(gdim, 3, 1, (ReferenceCellType::Triangle, 1));
builder.add_point(0, &[0.0, 0.0, 0.0]);
builder.add_point(1, &[1.0, 0.0, 0.0]);
builder.add_point(2, &[0.0, 1.0, 0.0]);
builder.add_cell(0, &[0, 1, 2]);

let grid = builder.create_grid();

Implementations

impl<T: Scalar> SingleElementGridBuilder<T>

pub fn new(gdim: usize, data: (ReferenceCellType, usize)) -> Self

Create a new grid builder

Examples found in repository

ndgrid/examples/test_parallel_io.rs (line 36)

fn example_single_element_grid<C: Communicator>(
    comm: &C,
    n: usize,
) -> ParallelGridImpl<'_, C, SingleElementGrid<f64, CiarletElement<f64, IdentityMap, f64>>> {
    let rank = comm.rank();

    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));

    if rank == 0 {
        create_single_element_grid_data(&mut b, n);
        b.create_parallel_grid_root(comm, GraphPartitioner::None)
    } else {
        b.create_parallel_grid(comm, 0)
    }
}

ndgrid/examples/test_parallel_grid.rs (line 12)

fn test_noncontiguous_numbering<C: Communicator>(comm: &C) {
    let rank = comm.rank();
    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));

    let g = if rank == 0 {
        let n = 5;
        for y in 0..n {
            for x in 0..n {
                b.add_point(
                    2 * (y * n + x) + 5,
                    &[x as f64 / (n - 1) as f64, y as f64 / (n - 1) as f64, 0.0],
                );
            }
        }

        for i in 0..n - 1 {
            for j in 0..n - 1 {
                b.add_cell(
                    3 * (i * (n - 1) + j),
                    &[
                        2 * (j * n + i) + 5,
                        2 * (j * n + i + 1) + 5,
                        2 * (j * n + i + n) + 5,
                        2 * (j * n + i + n + 1) + 5,
                    ],
                );
            }
        }

        b.create_parallel_grid_root(comm, GraphPartitioner::None)
    } else {
        b.create_parallel_grid(comm, 0)
    };

    assert!(g.local_grid().entity_count(ReferenceCellType::Point) > 0);
}

ndgrid/examples/single_element_grid.rs (line 13)

fn main() {
    // When creating the grid builder, we give the physical/geometric dimension (3) and the cell type
    // and degree of the element
    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
    // Add six points with ids 0 to 5
    b.add_point(0, &[0.0, 0.0, 0.0]);
    b.add_point(1, &[1.0, 0.0, 0.0]);
    b.add_point(2, &[2.0, 0.0, 0.2]);
    b.add_point(3, &[0.0, 1.0, 0.0]);
    b.add_point(4, &[1.0, 1.0, -0.2]);
    b.add_point(5, &[2.0, 1.0, 0.0]);
    // Add two cells
    b.add_cell(0, &[0, 1, 3, 4]);
    b.add_cell(1, &[1, 2, 4, 5]);
    // Create the grid
    let grid = b.create_grid();

    // Print the coordinates or each point in the mesh
    let mut coords = vec![0.0; grid.geometry_dim()];
    for point in grid.entity_iter(ReferenceCellType::Point) {
        point.geometry().points().collect::<Vec<_>>()[0].coords(coords.as_mut_slice());
        println!("point {}: {:#?}", point.local_index(), coords);
    }

    // Print the vertices of each cell
    for cell in grid.entity_iter(ReferenceCellType::Quadrilateral) {
        println!(
            "cell {}: {:?} ",
            cell.local_index(),
            cell.topology()
                .sub_entity_iter(ReferenceCellType::Point)
                .collect::<Vec<_>>()
        );
    }
}

ndgrid/examples/parallel_io.rs (line 29)

fn main() {
    let universe: Universe = mpi::initialize().unwrap();
    let comm = universe.world();
    let rank = comm.rank();

    let g = if rank == 0 {
        // Create a grid using the shapes module: unit_cube_boundary will mesh the surface of a cube
        let serial_g = shapes::unit_cube_boundary::<f64>(4, 5, 4, ReferenceCellType::Triangle);

        // Distribute this grid across processes
        serial_g.distribute(&comm, GraphPartitioner::None)
    } else {
        let b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
        b.create_parallel_grid(&comm, 0)
    };

    // If the serde option is used, the raw grid data can be exported in RON format
    g.export_as_ron("_unit_cube_boundary_parallel.ron");

    // Wait for export to finish
    comm.barrier();

    // A grid can be re-imported from raw RON data. Note that it must be imported on the same number of processes as it was exported using
    let g2 = ParallelGridImpl::<'_, _, SingleElementGrid::<f64, CiarletElement<f64, IdentityMap, f64>>>::import_from_ron(&comm, "_unit_cube_boundary_parallel.ron");

    // Print the first 5 cells of each grid on process 0
    if rank == 0 {
        println!("The first 5 cells of the grids");
        for (cell, cell2) in izip!(
            g.local_grid().entity_iter(ReferenceCellType::Triangle),
            g2.local_grid().entity_iter(ReferenceCellType::Triangle)
        )
        .take(5)
        {
            println!(
                "{:?} {:?}",
                cell.topology()
                    .sub_entity_iter(ReferenceCellType::Point)
                    .collect::<Vec<_>>(),
                cell2
                    .topology()
                    .sub_entity_iter(ReferenceCellType::Point)
                    .collect::<Vec<_>>(),
            );
        }
    }
}

ndgrid/examples/io.rs (line 53)

fn main() {
    // Create a grid using the shapes module: unit_cube_boundary will mesh the surface of a cube
    let g = shapes::unit_cube_boundary::<f64>(4, 5, 4, ReferenceCellType::Triangle);

    // If the serde option is used, the raw grid data can be exported in RON format
    g.export_as_ron("_unit_cube_boundary.ron");

    // A grid can be re-imported from raw RON data
    let g2 = SingleElementGrid::<f64, CiarletElement<f64, IdentityMap, f64>>::import_from_ron(
        "_unit_cube_boundary.ron",
    );

    // Print the first 5 cells of each grid
    println!("The first 5 cells of the grids");
    for (cell, cell2) in izip!(
        g.entity_iter(ReferenceCellType::Triangle),
        g2.entity_iter(ReferenceCellType::Triangle)
    )
    .take(5)
    {
        println!(
            "{:?} {:?}",
            cell.topology()
                .sub_entity_iter(ReferenceCellType::Point)
                .collect::<Vec<_>>(),
            cell2
                .topology()
                .sub_entity_iter(ReferenceCellType::Point)
                .collect::<Vec<_>>(),
        );
    }

    println!();

    // Alternatively, grids can be exported and imported to/from Gmsh files

    // Export the grid as a Gmsh .msh file
    g.export_as_gmsh("_unit_cube_boundary.msh");

    // To import from a Gmsh .msh file, a builder is used
    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
    b.import_from_gmsh("_unit_cube_boundary.msh");
    let g3 = b.create_grid();

    // Print the first 5 cells of each grid
    println!("The first 5 cells of the grids");
    for (cell, cell3) in izip!(
        g.entity_iter(ReferenceCellType::Triangle),
        g3.entity_iter(ReferenceCellType::Triangle)
    )
    .take(5)
    {
        println!(
            "{:?} {:?}",
            cell.topology()
                .sub_entity_iter(ReferenceCellType::Point)
                .collect::<Vec<_>>(),
            cell3
                .topology()
                .sub_entity_iter(ReferenceCellType::Point)
                .collect::<Vec<_>>(),
        );
    }
}

ndgrid/examples/test_push_forward.rs (line 16)

fn test_lagrange_push_forward() {
    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
    b.add_point(0, &[0.0, 0.0, 0.0]);
    b.add_point(1, &[1.0, 0.0, 0.0]);
    b.add_point(2, &[2.0, 0.0, 1.0]);
    b.add_point(3, &[0.0, 1.0, 0.0]);
    b.add_cell(0, &[0, 1, 3]);
    b.add_cell(1, &[1, 2, 3]);
    let grid = b.create_grid();

    let e = lagrange::create::<f64, f64>(ReferenceCellType::Triangle, 4, Continuity::Standard);

    let npts = 5;

    let mut cell0_points = rlst_dynamic_array!(f64, [2, npts]);
    for i in 0..npts {
        cell0_points[[1, i]] = i as f64 / (npts - 1) as f64;
        cell0_points[[0, i]] = 1.0 - cell0_points[[1, i]];
    }
    let mut cell0_table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
    e.tabulate(&cell0_points, 0, &mut cell0_table);

    let mut points = rlst_dynamic_array!(f64, [2, npts]);
    for i in 0..npts {
        points[[0, i]] = 0.0;
        points[[1, i]] = i as f64 / (npts - 1) as f64;
    }

    let mut table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
    e.tabulate(&points, 0, &mut table);

    let gmap = grid.geometry_map(ReferenceCellType::Triangle, 1, &points);

    let mut jacobians = rlst_dynamic_array!(f64, [grid.geometry_dim(), grid.topology_dim(), npts]);
    let mut jinv = rlst_dynamic_array!(f64, [grid.topology_dim(), grid.geometry_dim(), npts]);
    let mut jdets = vec![0.0; npts];

    gmap.jacobians_inverses_dets(1, &mut jacobians, &mut jinv, &mut jdets);

    let mut cell1_table = DynArray::<f64, 4>::from_shape(table.shape());
    e.push_forward(&table, 0, &jacobians, &jdets, &jinv, &mut cell1_table);

    // Check that basis functions are continuous between cells
    for (cell0_dof, cell1_dof) in izip!(
        e.entity_closure_dofs(1, 1).unwrap(),
        e.entity_closure_dofs(1, 0).unwrap()
    ) {
        for i in 0..npts {
            assert_relative_eq!(
                cell0_table[[0, i, *cell0_dof, 0]],
                cell1_table[[0, i, *cell1_dof, 0]],
                epsilon = 1e-10
            );
        }
    }
}

/// Test Ravaiart-Thomas push forward
fn test_rt_push_forward() {
    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
    b.add_point(0, &[0.0, 0.0, 0.0]);
    b.add_point(1, &[1.0, 0.0, 0.0]);
    b.add_point(2, &[2.0, 0.0, 1.0]);
    b.add_point(3, &[0.0, 1.0, 0.0]);
    b.add_cell(0, &[0, 1, 3]);
    b.add_cell(1, &[1, 2, 3]);
    let grid = b.create_grid();

    let e =
        raviart_thomas::create::<f64, f64>(ReferenceCellType::Triangle, 4, Continuity::Standard);

    let npts = 5;

    let mut cell0_points = rlst_dynamic_array!(f64, [2, npts]);
    for i in 0..npts {
        cell0_points[[1, i]] = i as f64 / (npts - 1) as f64;
        cell0_points[[0, i]] = 1.0 - cell0_points[[1, i]];
    }
    let mut cell0_table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
    e.tabulate(&cell0_points, 0, &mut cell0_table);

    let mut points = rlst_dynamic_array!(f64, [2, npts]);
    for i in 0..npts {
        points[[0, i]] = 0.0;
        points[[1, i]] = i as f64 / (npts - 1) as f64;
    }

    let mut table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
    e.tabulate(&points, 0, &mut table);

    let gmap = grid.geometry_map(ReferenceCellType::Triangle, 1, &points);

    let mut jacobians = rlst_dynamic_array!(f64, [grid.geometry_dim(), grid.topology_dim(), npts]);
    let mut jinv = rlst_dynamic_array!(f64, [grid.topology_dim(), grid.geometry_dim(), npts]);
    let mut jdets = vec![0.0; npts];

    gmap.jacobians_inverses_dets(1, &mut jacobians, &mut jinv, &mut jdets);

    let mut cell1_table =
        DynArray::<f64, 4>::from_shape([table.shape()[0], table.shape()[1], table.shape()[2], 3]);
    e.push_forward(&table, 0, &jacobians, &jdets, &jinv, &mut cell1_table);

    // Check that basis functions dotted with normal to edge are continuous between cells
    for (cell0_dof, cell1_dof) in izip!(
        e.entity_closure_dofs(1, 0).unwrap(),
        e.entity_closure_dofs(1, 1).unwrap()
    ) {
        for i in 0..npts {
            assert_relative_eq!(
                (cell0_table[[0, i, *cell0_dof, 0]] + cell0_table[[0, i, *cell0_dof, 1]])
                    / f64::sqrt(2.0),
                (cell1_table[[0, i, *cell1_dof, 0]]
                    + cell1_table[[0, i, *cell1_dof, 1]]
                    + 2.0 * cell1_table[[0, i, *cell1_dof, 2]])
                    / f64::sqrt(6.0),
                epsilon = 1e-10
            );
        }
    }
}

/// Test Nedelec push forward
fn test_nc_push_forward() {
    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
    b.add_point(0, &[0.0, 0.0, 0.0]);
    b.add_point(1, &[1.0, 0.0, 0.0]);
    b.add_point(2, &[2.0, 0.0, 1.0]);
    b.add_point(3, &[0.0, 1.0, 0.0]);
    b.add_cell(0, &[0, 1, 3]);
    b.add_cell(1, &[1, 2, 3]);
    let grid = b.create_grid();

    let e = nedelec::create::<f64, f64>(ReferenceCellType::Triangle, 4, Continuity::Standard);

    let npts = 5;

    let mut cell0_points = rlst_dynamic_array!(f64, [2, npts]);
    for i in 0..npts {
        cell0_points[[1, i]] = i as f64 / (npts - 1) as f64;
        cell0_points[[0, i]] = 1.0 - cell0_points[[1, i]];
    }
    let mut cell0_table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
    e.tabulate(&cell0_points, 0, &mut cell0_table);

    let mut points = rlst_dynamic_array!(f64, [2, npts]);
    for i in 0..npts {
        points[[0, i]] = 0.0;
        points[[1, i]] = i as f64 / (npts - 1) as f64;
    }

    let mut table = DynArray::<f64, 4>::from_shape(e.tabulate_array_shape(0, npts));
    e.tabulate(&points, 0, &mut table);

    let gmap = grid.geometry_map(ReferenceCellType::Triangle, 1, &points);

    let mut jacobians = rlst_dynamic_array!(f64, [grid.geometry_dim(), grid.topology_dim(), npts]);
    let mut jinv = rlst_dynamic_array!(f64, [grid.topology_dim(), grid.geometry_dim(), npts]);
    let mut jdets = vec![0.0; npts];

    gmap.jacobians_inverses_dets(1, &mut jacobians, &mut jinv, &mut jdets);

    let mut cell1_table =
        DynArray::<f64, 4>::from_shape([table.shape()[0], table.shape()[1], table.shape()[2], 3]);
    e.push_forward(&table, 0, &jacobians, &jdets, &jinv, &mut cell1_table);

    // Check that basis functions dotted with tangent to edge are continuous between cells
    for (cell0_dof, cell1_dof) in izip!(
        e.entity_closure_dofs(1, 0).unwrap(),
        e.entity_closure_dofs(1, 1).unwrap()
    ) {
        for i in 0..npts {
            assert_relative_eq!(
                (cell0_table[[0, i, *cell0_dof, 0]] - cell0_table[[0, i, *cell0_dof, 1]])
                    / f64::sqrt(2.0),
                (cell1_table[[0, i, *cell1_dof, 0]] - cell1_table[[0, i, *cell1_dof, 1]])
                    / f64::sqrt(2.0),
                epsilon = 1e-10
            );
        }
    }
}

Additional examples can be found in:

• ndgrid/examples/test_partitioners.rs
• ndgrid/examples/parallel_grid.rs

pub fn new_with_capacity( gdim: usize, npoints: usize, ncells: usize, data: (ReferenceCellType, usize), ) -> Self

Create a new grid builder with capacaty for a given number of points and cells

Trait Implementations

impl<T: Scalar> Builder for SingleElementGridBuilder<T>

type Grid = SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>

The type of the grid that the builder creates

type T = T

The floating point type used for coordinates

type CellData<'a> = &'a [usize]

The type of the data that is input to add a cell

type EntityDescriptor = ReferenceCellType

Type used as identifier of different entity types

fn add_point(&mut self, id: usize, data: &[T])

Add a point to the grid

fn add_cell(&mut self, id: usize, cell_data: &[usize])

Add a cell to the grid

fn add_cell_from_nodes_and_type( &mut self, id: usize, nodes: &[usize], cell_type: ReferenceCellType, cell_degree: usize, )

Add a cell to the grid

fn create_grid(&self) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>

Create the grid

fn point_count(&self) -> usize

Number of points

fn cell_count(&self) -> usize

Number of cells

fn point_indices_to_ids(&self) -> &[usize]

Get the insertion ids of each point

fn cell_indices_to_ids(&self) -> &[usize]

Get the insertion ids of each cell

fn cell_points(&self, index: usize) -> &[usize]

Get the indices of the points of a cell

fn cell_vertices(&self, index: usize) -> &[usize]

Get the indices of the points of a cell

fn point(&self, index: usize) -> &[T]

Get the coordinates of a point

fn points(&self) -> &[T]

Get all points

fn cell_type(&self, _index: usize) -> ReferenceCellType

Get the type of a cell

fn cell_degree(&self, _index: usize) -> usize

Get the degree of a cell’s geometry

fn gdim(&self) -> usize

Geometric dimension

fn tdim(&self) -> usize

Topoligical dimension

fn npts(&self, _cell_type: Self::EntityDescriptor, _degree: usize) -> usize

Number of points in a cell with the given type and degree

impl<T: Debug + Scalar> Debug for SingleElementGridBuilder<T>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.