xref: /libCEED/rust/libceed/src/lib.rs (revision bf55b007047d3d1ea1227af6b83a1b17e380b530)

// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors.
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
//
// SPDX-License-Identifier: BSD-2-Clause
//
// This file is part of CEED:  http://github.com/ceed

// Fenced `rust` code blocks included from README.md are executed as part of doctests.
#![doc = include_str!("../README.md")]
// -----------------------------------------------------------------------------
// Exceptions
// -----------------------------------------------------------------------------
#![allow(non_snake_case)]

// -----------------------------------------------------------------------------
// Crate prelude
// -----------------------------------------------------------------------------
use crate::prelude::*;
use std::sync::Once;

pub mod prelude {
    pub(crate) use libceed_sys::bind_ceed;
    pub(crate) use std::convert::TryFrom;
    pub(crate) use std::ffi::{CStr, CString};
    pub(crate) use std::fmt;
    pub(crate) use std::marker::PhantomData;
}

// -----------------------------------------------------------------------------
// Modules
// -----------------------------------------------------------------------------
pub mod basis;
pub mod elem_restriction;
pub mod operator;
pub mod qfunction;
pub mod vector;

// -----------------------------------------------------------------------------
// Typedef for scalar
// -----------------------------------------------------------------------------
pub type Scalar = bind_ceed::CeedScalar;

// -----------------------------------------------------------------------------
// Constants for library
// -----------------------------------------------------------------------------
const MAX_BUFFER_LENGTH: usize = 4096;
pub const MAX_QFUNCTION_FIELDS: usize = 16;
pub const CEED_STRIDES_BACKEND: [i32; 3] = [0; 3];
pub const EPSILON: crate::Scalar = bind_ceed::CEED_EPSILON as crate::Scalar;

// -----------------------------------------------------------------------------
// Enums for libCEED
// -----------------------------------------------------------------------------
#[derive(Clone, Copy, PartialEq, Eq)]
/// Many Ceed interfaces take or return pointers to memory.  This enum is used to
/// specify where the memory being provided or requested must reside.
pub enum MemType {
    Host = bind_ceed::CeedMemType_CEED_MEM_HOST as isize,
    Device = bind_ceed::CeedMemType_CEED_MEM_DEVICE as isize,
}

#[derive(Clone, Copy, PartialEq, Eq)]
// OwnPointer will not be used by user but is included for internal use
#[allow(dead_code)]
/// Conveys ownership status of arrays passed to Ceed interfaces.
pub(crate) enum CopyMode {
    CopyValues = bind_ceed::CeedCopyMode_CEED_COPY_VALUES as isize,
    UsePointer = bind_ceed::CeedCopyMode_CEED_USE_POINTER as isize,
    OwnPointer = bind_ceed::CeedCopyMode_CEED_OWN_POINTER as isize,
}

#[derive(Clone, Copy, PartialEq, Eq)]
/// Denotes type of vector norm to be computed
pub enum NormType {
    One = bind_ceed::CeedNormType_CEED_NORM_1 as isize,
    Two = bind_ceed::CeedNormType_CEED_NORM_2 as isize,
    Max = bind_ceed::CeedNormType_CEED_NORM_MAX as isize,
}

#[derive(Clone, Copy, PartialEq, Eq)]
/// Denotes whether a linear transformation or its transpose should be applied
pub enum TransposeMode {
    NoTranspose = bind_ceed::CeedTransposeMode_CEED_NOTRANSPOSE as isize,
    Transpose = bind_ceed::CeedTransposeMode_CEED_TRANSPOSE as isize,
}

#[derive(Clone, Copy, PartialEq, Eq)]
/// Type of quadrature; also used for location of nodes
pub enum QuadMode {
    Gauss = bind_ceed::CeedQuadMode_CEED_GAUSS as isize,
    GaussLobatto = bind_ceed::CeedQuadMode_CEED_GAUSS_LOBATTO as isize,
}

#[derive(Clone, Copy, PartialEq, Eq)]
/// Type of basis shape to create non-tensor H1 element basis
pub enum ElemTopology {
    Line = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_LINE as isize,
    Triangle = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_TRIANGLE as isize,
    Quad = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_QUAD as isize,
    Tet = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_TET as isize,
    Pyramid = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_PYRAMID as isize,
    Prism = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_PRISM as isize,
    Hex = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_HEX as isize,
}

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
/// Basis evaluation mode
pub enum EvalMode {
    None = bind_ceed::CeedEvalMode_CEED_EVAL_NONE as isize,
    Interp = bind_ceed::CeedEvalMode_CEED_EVAL_INTERP as isize,
    Grad = bind_ceed::CeedEvalMode_CEED_EVAL_GRAD as isize,
    Div = bind_ceed::CeedEvalMode_CEED_EVAL_DIV as isize,
    Curl = bind_ceed::CeedEvalMode_CEED_EVAL_CURL as isize,
    Weight = bind_ceed::CeedEvalMode_CEED_EVAL_WEIGHT as isize,
}
impl EvalMode {
    pub(crate) fn from_u32(value: u32) -> EvalMode {
        match value {
            bind_ceed::CeedEvalMode_CEED_EVAL_NONE => EvalMode::None,
            bind_ceed::CeedEvalMode_CEED_EVAL_INTERP => EvalMode::Interp,
            bind_ceed::CeedEvalMode_CEED_EVAL_GRAD => EvalMode::Grad,
            bind_ceed::CeedEvalMode_CEED_EVAL_DIV => EvalMode::Div,
            bind_ceed::CeedEvalMode_CEED_EVAL_CURL => EvalMode::Curl,
            bind_ceed::CeedEvalMode_CEED_EVAL_WEIGHT => EvalMode::Weight,
            _ => panic!("Unknown value: {}", value),
        }
    }
}

// -----------------------------------------------------------------------------
// Ceed error
// -----------------------------------------------------------------------------
pub type Result<T> = std::result::Result<T, Error>;

/// libCEED error messages - returning an Error without painc!ing indicates
///   the function call failed but the data is not corrupted
#[derive(Debug)]
pub struct Error {
    pub message: String,
}

impl fmt::Display for Error {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "{}", self.message)
    }
}

// -----------------------------------------------------------------------------
// Internal crate contents
// -----------------------------------------------------------------------------
pub use crate::{
    basis::{Basis, BasisOpt},
    elem_restriction::{ElemRestriction, ElemRestrictionOpt},
    operator::{CompositeOperator, Operator, OperatorField},
    qfunction::{
        QFunction, QFunctionByName, QFunctionField, QFunctionInputs, QFunctionOpt, QFunctionOutputs,
    },
    vector::{Vector, VectorOpt, VectorSliceWrapper},
};

// -----------------------------------------------------------------------------
// Internal error checker
// -----------------------------------------------------------------------------
#[doc(hidden)]
pub(crate) fn check_error<F>(ceed_ptr: F, ierr: i32) -> Result<i32>
where
    F: FnOnce() -> bind_ceed::Ceed,
{
    // Return early if code is clean
    if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS {
        return Ok(ierr);
    }
    // Retrieve error message
    let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut();
    let c_str = unsafe {
        bind_ceed::CeedGetErrorMessage(ceed_ptr(), &mut ptr);
        std::ffi::CStr::from_ptr(ptr)
    };
    let message = c_str.to_string_lossy().to_string();
    // Panic if negative code, otherwise return error
    if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS {
        panic!("{}", message);
    }
    Err(Error { message })
}

// -----------------------------------------------------------------------------
// Ceed error handler
// -----------------------------------------------------------------------------
pub enum ErrorHandler {
    ErrorAbort,
    ErrorExit,
    ErrorReturn,
    ErrorStore,
}

// -----------------------------------------------------------------------------
// Ceed context wrapper
// -----------------------------------------------------------------------------
/// A Ceed is a library context representing control of a logical hardware
/// resource.
#[derive(Debug)]
pub struct Ceed {
    ptr: bind_ceed::Ceed,
}

// -----------------------------------------------------------------------------
// Destructor
// -----------------------------------------------------------------------------
impl Drop for Ceed {
    fn drop(&mut self) {
        unsafe {
            bind_ceed::CeedDestroy(&mut self.ptr);
        }
    }
}

// -----------------------------------------------------------------------------
// Cloning
// -----------------------------------------------------------------------------
impl Clone for Ceed {
    /// Perform a shallow clone of a Ceed context
    ///
    /// ```
    /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial");
    /// let ceed_clone = ceed.clone();
    ///
    /// println!(" original:{} \n clone:{}", ceed, ceed_clone);
    /// ```
    fn clone(&self) -> Self {
        let mut ptr_clone = std::ptr::null_mut();
        self.check_error(unsafe { bind_ceed::CeedReferenceCopy(self.ptr, &mut ptr_clone) })
            .expect("failed to clone Ceed");
        Self { ptr: ptr_clone }
    }
}

// -----------------------------------------------------------------------------
// Display
// -----------------------------------------------------------------------------
impl fmt::Display for Ceed {
    /// View a Ceed
    ///
    /// ```
    /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial");
    /// println!("{}", ceed);
    /// ```
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        let mut ptr = std::ptr::null_mut();
        let mut sizeloc = crate::MAX_BUFFER_LENGTH;
        let cstring = unsafe {
            let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc);
            bind_ceed::CeedView(self.ptr, file);
            bind_ceed::fclose(file);
            CString::from_raw(ptr)
        };
        cstring.to_string_lossy().fmt(f)
    }
}

static REGISTER: Once = Once::new();

// -----------------------------------------------------------------------------
// Object constructors
// -----------------------------------------------------------------------------
impl Ceed {
    #[cfg_attr(feature = "katexit", katexit::katexit)]
    /// Returns a Ceed context initialized with the specified resource
    ///
    /// # arguments
    ///
    /// * `resource` - Resource to use, e.g., "/cpu/self"
    ///
    /// ```
    /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial");
    /// ```
    pub fn init(resource: &str) -> Self {
        Ceed::init_with_error_handler(resource, ErrorHandler::ErrorStore)
    }

    /// Returns a Ceed context initialized with the specified resource
    ///
    /// # arguments
    ///
    /// * `resource` - Resource to use, e.g., "/cpu/self"
    ///
    /// ```
    /// let ceed = libceed::Ceed::init_with_error_handler(
    ///     "/cpu/self/ref/serial",
    ///     libceed::ErrorHandler::ErrorAbort,
    /// );
    /// ```
    pub fn init_with_error_handler(resource: &str, handler: ErrorHandler) -> Self {
        REGISTER.call_once(|| unsafe {
            bind_ceed::CeedRegisterAll();
            bind_ceed::CeedQFunctionRegisterAll();
        });

        // Convert to C string
        let c_resource = CString::new(resource).expect("CString::new failed");

        // Get error handler pointer
        let eh = match handler {
            ErrorHandler::ErrorAbort => bind_ceed::CeedErrorAbort,
            ErrorHandler::ErrorExit => bind_ceed::CeedErrorExit,
            ErrorHandler::ErrorReturn => bind_ceed::CeedErrorReturn,
            ErrorHandler::ErrorStore => bind_ceed::CeedErrorStore,
        };

        // Call to libCEED
        let mut ptr = std::ptr::null_mut();
        let mut ierr = unsafe { bind_ceed::CeedInit(c_resource.as_ptr(), &mut ptr) };
        if ierr != 0 {
            panic!("Error initializing backend resource: {}", resource)
        }
        ierr = unsafe { bind_ceed::CeedSetErrorHandler(ptr, Some(eh)) };
        let ceed = Ceed { ptr };
        ceed.check_error(ierr).unwrap();
        ceed
    }

    /// Default initializer for testing
    #[doc(hidden)]
    pub fn default_init() -> Self {
        // Convert to C string
        let resource = "/cpu/self/ref/serial";
        crate::Ceed::init(resource)
    }

    /// Internal error checker
    #[doc(hidden)]
    fn check_error(&self, ierr: i32) -> Result<i32> {
        // Return early if code is clean
        if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS {
            return Ok(ierr);
        }
        // Retrieve error message
        let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut();
        let c_str = unsafe {
            bind_ceed::CeedGetErrorMessage(self.ptr, &mut ptr);
            std::ffi::CStr::from_ptr(ptr)
        };
        let message = c_str.to_string_lossy().to_string();
        // Panic if negative code, otherwise return error
        if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS {
            panic!("{}", message);
        }
        Err(Error { message })
    }

    /// Returns full resource name for a Ceed context
    ///
    /// ```
    /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial");
    /// let resource = ceed.resource();
    ///
    /// assert_eq!(resource, "/cpu/self/ref/serial".to_string())
    /// ```
    pub fn resource(&self) -> String {
        let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut();
        let c_str = unsafe {
            bind_ceed::CeedGetResource(self.ptr, &mut ptr);
            std::ffi::CStr::from_ptr(ptr)
        };
        c_str.to_string_lossy().to_string()
    }

    /// Returns a Vector of the specified length (does not allocate memory)
    ///
    /// # arguments
    ///
    /// * `n` - Length of vector
    ///
    /// ```
    /// # use libceed::prelude::*;
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let vec = ceed.vector(10)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn vector<'a>(&self, n: usize) -> Result<Vector<'a>> {
        Vector::create(self, n)
    }

    /// Create a Vector initialized with the data (copied) from a slice
    ///
    /// # arguments
    ///
    /// * `slice` - Slice containing data
    ///
    /// ```
    /// # use libceed::prelude::*;
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let vec = ceed.vector_from_slice(&[1., 2., 3.])?;
    /// assert_eq!(vec.length(), 3);
    /// # Ok(())
    /// # }
    /// ```
    pub fn vector_from_slice<'a>(&self, slice: &[crate::Scalar]) -> Result<Vector<'a>> {
        Vector::from_slice(self, slice)
    }

    /// Returns an ElemRestriction, $\mathcal{E}$, which extracts the degrees of
    ///   freedom for each element from the local vector into the element vector
    ///   or assembles contributions from each element in the element vector to
    ///   the local vector
    ///
    /// # arguments
    ///
    /// * `nelem`      - Number of elements described in the offsets array
    /// * `elemsize`   - Size (number of "nodes") per element
    /// * `ncomp`      - Number of field components per interpolation node (1
    ///                    for scalar fields)
    /// * `compstride` - Stride between components for the same Lvector "node".
    ///                    Data for node `i`, component `j`, element `k` can be
    ///                    found in the Lvector at index
    ///                    `offsets[i + k*elemsize] + j*compstride`.
    /// * `lsize`      - The size of the Lvector. This vector may be larger
    ///                    than the elements and fields given by this
    ///                    restriction.
    /// * `mtype`      - Memory type of the offsets array, see CeedMemType
    /// * `offsets`    - Array of shape `[nelem, elemsize]`. Row `i` holds the
    ///                    ordered list of the offsets (into the input CeedVector)
    ///                    for the unknowns corresponding to element `i`, where
    ///                    `0 <= i < nelem`. All offsets must be in the range
    ///                    `[0, lsize - 1]`.
    ///
    /// ```
    /// # use libceed::{prelude::*, MemType};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let nelem = 3;
    /// let mut ind: Vec<i32> = vec![0; 2 * nelem];
    /// for i in 0..nelem {
    ///     ind[2 * i + 0] = i as i32;
    ///     ind[2 * i + 1] = (i + 1) as i32;
    /// }
    /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn elem_restriction<'a>(
        &self,
        nelem: usize,
        elemsize: usize,
        ncomp: usize,
        compstride: usize,
        lsize: usize,
        mtype: MemType,
        offsets: &[i32],
    ) -> Result<ElemRestriction<'a>> {
        ElemRestriction::create(
            self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets,
        )
    }

    /// Returns an oriented ElemRestriction, $\mathcal{E}$, which extracts the
    ///   degrees of freedom for each element from the local vector into the
    ///   element vector or assembles contributions from each element in the
    ///   element vector to the local vector
    ///
    /// # arguments
    ///
    /// * `nelem`      - Number of elements described in the offsets array
    /// * `elemsize`   - Size (number of "nodes") per element
    /// * `ncomp`      - Number of field components per interpolation node (1
    ///                    for scalar fields)
    /// * `compstride` - Stride between components for the same Lvector "node".
    ///                    Data for node `i`, component `j`, element `k` can be
    ///                    found in the Lvector at index
    ///                    `offsets[i + k*elemsize] + j*compstride`.
    /// * `lsize`      - The size of the Lvector. This vector may be larger
    ///                    than the elements and fields given by this
    ///                    restriction.
    /// * `mtype`      - Memory type of the offsets array, see CeedMemType
    /// * `offsets`    - Array of shape `[nelem, elemsize]`. Row `i` holds the
    ///                    ordered list of the offsets (into the input CeedVector)
    ///                    for the unknowns corresponding to element `i`, where
    ///                    `0 <= i < nelem`. All offsets must be in the range
    ///                    `[0, lsize - 1]`.
    /// * `orients`    - Array of shape `[nelem, elemsize]`. Row `i` holds the
    ///                    ordered list of the orientations for the unknowns
    ///                    corresponding to element `i`, with bool `false` used
    ///                    for positively oriented and `true` to flip the
    ///                    orientation.
    ///
    /// ```
    /// # use libceed::{prelude::*, MemType};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let nelem = 3;
    /// let mut ind: Vec<i32> = vec![0; 2 * nelem];
    /// let mut orients: Vec<bool> = vec![false; 2 * nelem];
    /// for i in 0..nelem {
    ///     ind[2 * i + 0] = i as i32;
    ///     ind[2 * i + 1] = (i + 1) as i32;
    ///     orients[2 * i + 0] = (i % 2) > 0; // flip the dofs on element 1, 3, ...
    ///     orients[2 * i + 1] = (i % 2) > 0;
    /// }
    /// let r =
    ///     ceed.oriented_elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind, &orients)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn oriented_elem_restriction<'a>(
        &self,
        nelem: usize,
        elemsize: usize,
        ncomp: usize,
        compstride: usize,
        lsize: usize,
        mtype: MemType,
        offsets: &[i32],
        orients: &[bool],
    ) -> Result<ElemRestriction<'a>> {
        ElemRestriction::create_oriented(
            self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, orients,
        )
    }

    /// Returns a curl-oriented ElemRestriction, $\mathcal{E}$, which extracts
    ///   the degrees of freedom for each element from the local vector into
    ///   the element vector or assembles contributions from each element in
    ///   the element vector to the local vector
    ///
    /// # arguments
    ///
    /// * `nelem`       - Number of elements described in the offsets array
    /// * `elemsize`    - Size (number of "nodes") per element
    /// * `ncomp`       - Number of field components per interpolation node (1
    ///                     for scalar fields)
    /// * `compstride`  - Stride between components for the same Lvector "node".
    ///                     Data for node `i`, component `j`, element `k` can be
    ///                     found in the Lvector at index
    ///                     `offsets[i + k*elemsize] + j*compstride`.
    /// * `lsize`       - The size of the Lvector. This vector may be larger
    ///                     than the elements and fields given by this
    ///                     restriction.
    /// * `mtype`       - Memory type of the offsets array, see CeedMemType
    /// * `offsets`     - Array of shape `[nelem, elemsize]`. Row `i` holds the
    ///                     ordered list of the offsets (into the input CeedVector)
    ///                     for the unknowns corresponding to element `i`, where
    ///                     `0 <= i < nelem`. All offsets must be in the range
    ///                     `[0, lsize - 1]`.
    /// * `curlorients` - Array of shape `[nelem, 3 * elemsize]`. Row `i` holds
    ///                     a row-major tridiagonal matrix (`curlorients[i, 0] =
    ///                     curlorients[i, 3 * elemsize - 1] = 0`, where
    ///                     `0 <= i < nelem`) which is applied to the element
    ///                     unknowns upon restriction.
    ///
    /// ```
    /// # use libceed::{prelude::*, MemType};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let nelem = 3;
    /// let mut ind: Vec<i32> = vec![0; 2 * nelem];
    /// let mut curlorients: Vec<i8> = vec![0; 3 * 2 * nelem];
    /// for i in 0..nelem {
    ///     ind[2 * i + 0] = i as i32;
    ///     ind[2 * i + 1] = (i + 1) as i32;
    ///     curlorients[3 * 2 * i] = 0 as i8;
    ///     curlorients[3 * 2 * (i + 1) - 1] = 0 as i8;
    ///     if (i % 2 > 0) {
    ///         // T = [0  -1]
    ///         //     [-1  0]
    ///         curlorients[3 * 2 * i + 1] = 0 as i8;
    ///         curlorients[3 * 2 * i + 2] = -1 as i8;
    ///         curlorients[3 * 2 * i + 3] = -1 as i8;
    ///         curlorients[3 * 2 * i + 4] = 0 as i8;
    ///     } else {
    ///         // T = I
    ///         curlorients[3 * 2 * i + 1] = 1 as i8;
    ///         curlorients[3 * 2 * i + 2] = 0 as i8;
    ///         curlorients[3 * 2 * i + 3] = 0 as i8;
    ///         curlorients[3 * 2 * i + 4] = 1 as i8;
    ///     }
    /// }
    /// let r = ceed.curl_oriented_elem_restriction(
    ///     nelem,
    ///     2,
    ///     1,
    ///     1,
    ///     nelem + 1,
    ///     MemType::Host,
    ///     &ind,
    ///     &curlorients,
    /// )?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn curl_oriented_elem_restriction<'a>(
        &self,
        nelem: usize,
        elemsize: usize,
        ncomp: usize,
        compstride: usize,
        lsize: usize,
        mtype: MemType,
        offsets: &[i32],
        curlorients: &[i8],
    ) -> Result<ElemRestriction<'a>> {
        ElemRestriction::create_curl_oriented(
            self,
            nelem,
            elemsize,
            ncomp,
            compstride,
            lsize,
            mtype,
            offsets,
            curlorients,
        )
    }

    /// Returns an ElemRestriction, $\mathcal{E}$, from an local vector to
    ///   an element vector where data can be indexed from the `strides` array
    ///
    /// # arguments
    ///
    /// * `nelem`      - Number of elements described in the offsets array
    /// * `elemsize`   - Size (number of "nodes") per element
    /// * `ncomp`      - Number of field components per interpolation node (1
    ///                    for scalar fields)
    /// * `lsize`      - The size of the Lvector. This vector may be larger
    ///   than the elements and fields given by this restriction.
    /// * `strides`   - Array for strides between `[nodes, components, elements]`.
    ///                   Data for node `i`, component `j`, element `k` can be
    ///                   found in the Lvector at index
    ///                   `i*strides[0] + j*strides[1] + k*strides[2]`.
    ///                   CEED_STRIDES_BACKEND may be used with vectors created
    ///                   by a Ceed backend.
    ///
    /// ```
    /// # use libceed::prelude::*;
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let nelem = 3;
    /// let strides: [i32; 3] = [1, 2, 2];
    /// let r = ceed.strided_elem_restriction(nelem, 2, 1, nelem * 2, strides)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn strided_elem_restriction<'a>(
        &self,
        nelem: usize,
        elemsize: usize,
        ncomp: usize,
        lsize: usize,
        strides: [i32; 3],
    ) -> Result<ElemRestriction<'a>> {
        ElemRestriction::create_strided(self, nelem, elemsize, ncomp, lsize, strides)
    }

    /// Returns an $H^1$ tensor-product Basis
    ///
    /// # arguments
    ///
    /// * `dim`       - Topological dimension of element
    /// * `ncomp`     - Number of field components (1 for scalar fields)
    /// * `P1d`       - Number of Gauss-Lobatto nodes in one dimension.  The
    ///                   polynomial degree of the resulting `Q_k` element is
    ///                   `k=P-1`.
    /// * `Q1d`       - Number of quadrature points in one dimension
    /// * `interp1d`  - Row-major `(Q1d * P1d)` matrix expressing the values of
    ///                   nodal basis functions at quadrature points
    /// * `grad1d`    - Row-major `(Q1d * P1d)` matrix expressing derivatives of
    ///                   nodal basis functions at quadrature points
    /// * `qref1d`    - Array of length `Q1d` holding the locations of quadrature
    ///                   points on the 1D reference element `[-1, 1]`
    /// * `qweight1d` - Array of length `Q1d` holding the quadrature weights on
    ///                   the reference element
    ///
    /// ```
    /// # use libceed::prelude::*;
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let interp1d  = [ 0.62994317,  0.47255875, -0.14950343,  0.04700152,
    ///                  -0.07069480,  0.97297619,  0.13253993, -0.03482132,
    ///                  -0.03482132,  0.13253993,  0.97297619, -0.07069480,
    ///                   0.04700152, -0.14950343,  0.47255875,  0.62994317];
    /// let grad1d    = [-2.34183742,  2.78794489, -0.63510411,  0.18899664,
    ///                  -0.51670214, -0.48795249,  1.33790510, -0.33325047,
    //                    0.33325047, -1.33790510,  0.48795249,  0.51670214,
    ///                  -0.18899664,  0.63510411, -2.78794489,  2.34183742];
    /// let qref1d    = [-0.86113631, -0.33998104,  0.33998104,  0.86113631];
    /// let qweight1d = [ 0.34785485,  0.65214515,  0.65214515,  0.34785485];
    /// let b = ceed.
    /// basis_tensor_H1(2, 1, 4, 4, &interp1d, &grad1d, &qref1d, &qweight1d)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn basis_tensor_H1<'a>(
        &self,
        dim: usize,
        ncomp: usize,
        P1d: usize,
        Q1d: usize,
        interp1d: &[crate::Scalar],
        grad1d: &[crate::Scalar],
        qref1d: &[crate::Scalar],
        qweight1d: &[crate::Scalar],
    ) -> Result<Basis<'a>> {
        Basis::create_tensor_H1(
            self, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d,
        )
    }

    /// Returns an $H^1$ Lagrange tensor-product Basis
    ///
    /// # arguments
    ///
    /// * `dim`   - Topological dimension of element
    /// * `ncomp` - Number of field components (1 for scalar fields)
    /// * `P`     - Number of Gauss-Lobatto nodes in one dimension.  The
    ///               polynomial degree of the resulting `Q_k` element is `k=P-1`.
    /// * `Q`     - Number of quadrature points in one dimension
    /// * `qmode` - Distribution of the `Q` quadrature points (affects order of
    ///               accuracy for the quadrature)
    ///
    /// ```
    /// # use libceed::{prelude::*, QuadMode};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let b = ceed.basis_tensor_H1_Lagrange(2, 1, 3, 4, QuadMode::Gauss)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn basis_tensor_H1_Lagrange<'a>(
        &self,
        dim: usize,
        ncomp: usize,
        P: usize,
        Q: usize,
        qmode: QuadMode,
    ) -> Result<Basis<'a>> {
        Basis::create_tensor_H1_Lagrange(self, dim, ncomp, P, Q, qmode)
    }

    /// Returns an $H-1$ Basis
    ///
    /// # arguments
    ///
    /// * `topo`    - Topology of element, e.g. hypercube, simplex, ect
    /// * `ncomp`   - Number of field components (1 for scalar fields)
    /// * `nnodes`  - Total number of nodes
    /// * `nqpts`   - Total number of quadrature points
    /// * `interp`  - Row-major `(nqpts * nnodes)` matrix expressing the values of
    ///                 nodal basis functions at quadrature points
    /// * `grad`    - Row-major `(dim * nqpts * nnodes)` matrix expressing
    ///                 derivatives of nodal basis functions at quadrature points
    /// * `qref`    - Array of length `nqpts` holding the locations of quadrature
    ///                 points on the reference element `[-1, 1]`
    /// * `qweight` - Array of length `nqpts` holding the quadrature weights on
    ///                 the reference element
    ///
    /// ```
    /// # use libceed::{prelude::*, ElemTopology};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let interp = [
    ///     0.12000000,
    ///     0.48000000,
    ///     -0.12000000,
    ///     0.48000000,
    ///     0.16000000,
    ///     -0.12000000,
    ///     -0.12000000,
    ///     0.48000000,
    ///     0.12000000,
    ///     0.16000000,
    ///     0.48000000,
    ///     -0.12000000,
    ///     -0.11111111,
    ///     0.44444444,
    ///     -0.11111111,
    ///     0.44444444,
    ///     0.44444444,
    ///     -0.11111111,
    ///     -0.12000000,
    ///     0.16000000,
    ///     -0.12000000,
    ///     0.48000000,
    ///     0.48000000,
    ///     0.12000000,
    /// ];
    /// let grad = [
    ///     -1.40000000,
    ///     1.60000000,
    ///     -0.20000000,
    ///     -0.80000000,
    ///     0.80000000,
    ///     0.00000000,
    ///     0.20000000,
    ///     -1.60000000,
    ///     1.40000000,
    ///     -0.80000000,
    ///     0.80000000,
    ///     0.00000000,
    ///     -0.33333333,
    ///     0.00000000,
    ///     0.33333333,
    ///     -1.33333333,
    ///     1.33333333,
    ///     0.00000000,
    ///     0.20000000,
    ///     0.00000000,
    ///     -0.20000000,
    ///     -2.40000000,
    ///     2.40000000,
    ///     0.00000000,
    ///     -1.40000000,
    ///     -0.80000000,
    ///     0.00000000,
    ///     1.60000000,
    ///     0.80000000,
    ///     -0.20000000,
    ///     0.20000000,
    ///     -2.40000000,
    ///     0.00000000,
    ///     0.00000000,
    ///     2.40000000,
    ///     -0.20000000,
    ///     -0.33333333,
    ///     -1.33333333,
    ///     0.00000000,
    ///     0.00000000,
    ///     1.33333333,
    ///     0.33333333,
    ///     0.20000000,
    ///     -0.80000000,
    ///     0.00000000,
    ///     -1.60000000,
    ///     0.80000000,
    ///     1.40000000,
    /// ];
    /// let qref = [
    ///     0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333,
    ///     0.60000000,
    /// ];
    /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667];
    /// let b = ceed.basis_H1(
    ///     ElemTopology::Triangle,
    ///     1,
    ///     6,
    ///     4,
    ///     &interp,
    ///     &grad,
    ///     &qref,
    ///     &qweight,
    /// )?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn basis_H1<'a>(
        &self,
        topo: ElemTopology,
        ncomp: usize,
        nnodes: usize,
        nqpts: usize,
        interp: &[crate::Scalar],
        grad: &[crate::Scalar],
        qref: &[crate::Scalar],
        qweight: &[crate::Scalar],
    ) -> Result<Basis<'a>> {
        Basis::create_H1(
            self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight,
        )
    }

    /// Returns an $H(div)$ Basis
    ///
    /// # arguments
    ///
    /// * `topo`    - Topology of element, e.g. hypercube, simplex, ect
    /// * `ncomp`   - Number of field components (1 for scalar fields)
    /// * `nnodes`  - Total number of nodes
    /// * `nqpts`   - Total number of quadrature points
    /// * `interp`  - Row-major `(dim * nqpts * nnodes)` matrix expressing the
    ///                 values of basis functions at quadrature points
    /// * `div`     - Row-major `(nqpts * nnodes)` matrix expressing the
    ///                 divergence of basis functions at quadrature points
    /// * `qref`    - Array of length `nqpts` holding the locations of quadrature
    ///                 points on the reference element `[-1, 1]`
    /// * `qweight` - Array of length `nqpts` holding the quadrature weights on
    ///                 the reference element
    ///
    /// ```
    /// # use libceed::{prelude::*, ElemTopology};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let interp = [
    ///     0.00000000,
    ///     -1.57735027,
    ///     0.57735027,
    ///     0.00000000,
    ///     0.00000000,
    ///     -1.57735027,
    ///     0.57735027,
    ///     0.00000000,
    ///     0.00000000,
    ///     -1.57735027,
    ///     0.57735027,
    ///     0.00000000,
    ///     0.00000000,
    ///     -0.42264973,
    ///     -0.57735027,
    ///     0.00000000,
    ///     0.42264973,
    ///     0.00000000,
    ///     0.00000000,
    ///     0.57735027,
    ///     0.42264973,
    ///     0.00000000,
    ///     0.00000000,
    ///     0.57735027,
    ///     1.57735027,
    ///     0.00000000,
    ///     0.00000000,
    ///     -0.57735027,
    ///     1.57735027,
    ///     0.00000000,
    ///     0.00000000,
    ///     -0.57735027,
    /// ];
    /// let div = [
    ///     -1.00000000,
    ///     1.00000000,
    ///     -1.00000000,
    ///     1.00000000,
    ///     -1.00000000,
    ///     1.00000000,
    ///     -1.00000000,
    ///     1.00000000,
    ///     -1.00000000,
    ///     1.00000000,
    ///     -1.00000000,
    ///     1.00000000,
    ///     -1.00000000,
    ///     1.00000000,
    ///     -1.00000000,
    ///     1.00000000,
    /// ];
    /// let qref = [
    ///     0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026,
    ///     0.57735026,
    /// ];
    /// let qweight = [1.00000000, 1.00000000, 1.00000000, 1.00000000];
    /// let b = ceed.basis_Hdiv(ElemTopology::Quad, 1, 4, 4, &interp, &div, &qref, &qweight)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn basis_Hdiv<'a>(
        &self,
        topo: ElemTopology,
        ncomp: usize,
        nnodes: usize,
        nqpts: usize,
        interp: &[crate::Scalar],
        div: &[crate::Scalar],
        qref: &[crate::Scalar],
        qweight: &[crate::Scalar],
    ) -> Result<Basis<'a>> {
        Basis::create_Hdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight)
    }

    /// Returns an $H(curl)$ Basis
    ///
    /// # arguments
    ///
    /// * `topo`    - Topology of element, e.g. hypercube, simplex, ect
    /// * `ncomp`   - Number of field components (1 for scalar fields)
    /// * `nnodes`  - Total number of nodes
    /// * `nqpts`   - Total number of quadrature points
    /// * `interp`  - Row-major `(dim * nqpts * nnodes)` matrix expressing the
    ///                 values of basis functions at quadrature points
    /// * `curl`    - Row-major `(curl_comp * nqpts * nnodes)`, `curl_comp = 1 if
    ///                 dim < 3 else dim` matrix expressing the curl of basis
    ///                 functions at quadrature points
    /// * `qref`    - Array of length `nqpts` holding the locations of quadrature
    ///                 points on the reference element `[-1, 1]`
    /// * `qweight` - Array of length `nqpts` holding the quadrature weights on
    ///                 the reference element
    ///
    /// ```
    /// # use libceed::{prelude::*, ElemTopology};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let interp = [
    ///     -0.20000000,
    ///     0.20000000,
    ///     0.80000000,
    ///     -0.20000000,
    ///     0.20000000,
    ///     0.80000000,
    ///     -0.33333333,
    ///     0.33333333,
    ///     0.66666667,
    ///     -0.60000000,
    ///     0.60000000,
    ///     0.40000000,
    ///     0.20000000,
    ///     0.80000000,
    ///     0.20000000,
    ///     0.60000000,
    ///     0.40000000,
    ///     0.60000000,
    ///     0.33333333,
    ///     0.66666667,
    ///     0.33333333,
    ///     0.20000000,
    ///     0.80000000,
    ///     0.20000000,
    /// ];
    /// let curl = [
    ///     2.00000000,
    ///     -2.00000000,
    ///     -2.00000000,
    ///     2.00000000,
    ///     -2.00000000,
    ///     -2.00000000,
    ///     2.00000000,
    ///     -2.00000000,
    ///     -2.00000000,
    ///     2.00000000,
    ///     -2.00000000,
    ///     -2.00000000,
    /// ];
    /// let qref = [
    ///     0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333,
    ///     0.60000000,
    /// ];
    /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667];
    /// let b = ceed.basis_Hcurl(
    ///     ElemTopology::Triangle,
    ///     1,
    ///     3,
    ///     4,
    ///     &interp,
    ///     &curl,
    ///     &qref,
    ///     &qweight,
    /// )?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn basis_Hcurl<'a>(
        &self,
        topo: ElemTopology,
        ncomp: usize,
        nnodes: usize,
        nqpts: usize,
        interp: &[crate::Scalar],
        curl: &[crate::Scalar],
        qref: &[crate::Scalar],
        qweight: &[crate::Scalar],
    ) -> Result<Basis<'a>> {
        Basis::create_Hcurl(
            self, topo, ncomp, nnodes, nqpts, interp, curl, qref, qweight,
        )
    }

    /// Returns a QFunction for evaluating interior (volumetric) terms
    ///   of a weak form corresponding to the $L^2$ inner product
    ///
    /// $$
    /// \langle v, F(u) \rangle = \int_\Omega v \cdot f_0 \left( u, \nabla u \right) + \left( \nabla v \right) : f_1 \left( u, \nabla u \right),
    /// $$
    ///
    /// where $v \cdot f_0$ represents contraction over fields and $\nabla v : f_1$
    ///   represents contraction over both fields and spatial dimensions.
    ///
    /// # arguments
    ///
    /// * `vlength` - Vector length. Caller must ensure that number of
    ///                 quadrature points is a multiple of vlength.
    /// * `f`       - Boxed closure to evaluate weak form at quadrature points.
    ///
    /// ```
    /// # use libceed::{prelude::*, QFunctionInputs, QFunctionOutputs};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| {
    ///     // Iterate over quadrature points
    ///     v.iter_mut()
    ///         .zip(u.iter().zip(weights.iter()))
    ///         .for_each(|(v, (u, w))| *v = u * w);
    ///
    ///     // Return clean error code
    ///     0
    /// };
    ///
    /// let qf = ceed.q_function_interior(1, Box::new(user_f))?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn q_function_interior<'a>(
        &self,
        vlength: usize,
        f: Box<qfunction::QFunctionUserClosure>,
    ) -> Result<QFunction<'a>> {
        QFunction::create(self, vlength, f)
    }

    /// Returns a QFunction for evaluating interior (volumetric) terms
    ///   created by name
    ///
    /// # arguments
    ///
    /// * `name` - name of QFunction from libCEED gallery
    ///
    /// ```
    /// # use libceed::prelude::*;
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn q_function_interior_by_name<'a>(&self, name: &str) -> Result<QFunctionByName<'a>> {
        QFunctionByName::create(self, name)
    }

    /// Returns an Operator and associate a QFunction. A Basis and
    ///   ElemRestriction can be associated with QFunction fields via
    ///   set_field().
    ///
    /// # arguments
    ///
    /// * `qf`   - QFunction defining the action of the operator at quadrature
    ///              points
    /// * `dqf`  - QFunction defining the action of the Jacobian of the qf (or
    ///              qfunction_none)
    /// * `dqfT` - QFunction defining the action of the transpose of the
    ///              Jacobian of the qf (or qfunction_none)
    ///
    /// ```
    /// # use libceed::{prelude::*, QFunctionOpt};
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?;
    /// let op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn operator<'a, 'b>(
        &self,
        qf: impl Into<QFunctionOpt<'b>>,
        dqf: impl Into<QFunctionOpt<'b>>,
        dqfT: impl Into<QFunctionOpt<'b>>,
    ) -> Result<Operator<'a>> {
        Operator::create(self, qf, dqf, dqfT)
    }

    /// Returns an Operator that composes the action of several Operators
    ///
    /// ```
    /// # use libceed::prelude::*;
    /// # fn main() -> libceed::Result<()> {
    /// # let ceed = libceed::Ceed::default_init();
    /// let op = ceed.composite_operator()?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn composite_operator<'a>(&self) -> Result<CompositeOperator<'a>> {
        CompositeOperator::create(self)
    }
}

// -----------------------------------------------------------------------------
// Tests
// -----------------------------------------------------------------------------
#[cfg(test)]
mod tests {
    use super::*;

    fn ceed_t501() -> Result<()> {
        let resource = "/cpu/self/ref/blocked";
        let ceed = Ceed::init(resource);
        let nelem = 4;
        let p = 3;
        let q = 4;
        let ndofs = p * nelem - nelem + 1;

        // Vectors
        let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?;
        let mut qdata = ceed.vector(nelem * q)?;
        qdata.set_value(0.0)?;
        let mut u = ceed.vector(ndofs)?;
        u.set_value(1.0)?;
        let mut v = ceed.vector(ndofs)?;
        v.set_value(0.0)?;

        // Restrictions
        let mut indx: Vec<i32> = vec![0; 2 * nelem];
        for i in 0..nelem {
            for j in 0..2 {
                indx[2 * i + j] = (i + j) as i32;
            }
        }
        let rx = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &indx)?;
        let mut indu: Vec<i32> = vec![0; p * nelem];
        for i in 0..nelem {
            for j in 0..p {
                indu[p * i + j] = (i + j) as i32;
            }
        }
        let ru = ceed.elem_restriction(nelem, p, 1, 1, ndofs, MemType::Host, &indu)?;
        let strides: [i32; 3] = [1, q as i32, q as i32];
        let rq = ceed.strided_elem_restriction(nelem, q, 1, q * nelem, strides)?;

        // Bases
        let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?;
        let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?;

        // Build quadrature data
        let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?;
        ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)?
            .field("dx", &rx, &bx, VectorOpt::Active)?
            .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)?
            .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?
            .apply(&x, &mut qdata)?;

        // Mass operator
        let qf_mass = ceed.q_function_interior_by_name("MassApply")?;
        let op_mass = ceed
            .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)?
            .field("u", &ru, &bu, VectorOpt::Active)?
            .field("qdata", &rq, BasisOpt::None, &qdata)?
            .field("v", &ru, &bu, VectorOpt::Active)?
            .check()?;

        v.set_value(0.0)?;
        op_mass.apply(&u, &mut v)?;

        // Check
        let sum: Scalar = v.view()?.iter().sum();
        let error: Scalar = (sum - 2.0).abs();
        assert!(
            error < 50.0 * EPSILON,
            "Incorrect interval length computed. Expected: 2.0, Found: {}, Error: {:.12e}",
            sum,
            error
        );
        Ok(())
    }

    #[test]
    fn test_ceed_t501() {
        assert!(ceed_t501().is_ok());
    }
}

// -----------------------------------------------------------------------------