src/dolfinx_access.rs

/*!
 Utility functions accessing C++ side Dolfinx.
*/
use alg_tools::error::DynResult;
use alg_tools::fe_model::{
    base::RealLocalModel,
    p2_local_model::{P2LocalModel, Simplex},
};
use alg_tools::loc::Loc;
use anyhow::bail;
use pyo3::ffi::PyImport_AppendInittab;
use std::ffi::CString;
use std::mem::MaybeUninit;

/// Import helper required for the linking to the sys crates to actually happen.
#[allow(unused_imports)]
mod dummy_import {
    use dolfinx_sys;
    use nanobind_sys;
}

mod function;
pub use function::DolfinxPyFunction_f64;

/// C++ side routines
#[allow(dead_code)]
#[cxx::bridge(namespace = "dolfinx_access")]
pub mod ffi {
    /// Structure for coordinte-value pairs
    #[derive(Copy, Clone, Debug)]
    struct CoordValuePair {
        /// 2D coordinate
        x: [f64; 2],
        /// Value
        v: f64,
    }

    /// Structure for information about a dolfinx `Function`
    #[derive(Copy, Clone, Debug)]
    struct FunctionInfo {
        /// Domain dimension
        domain_dim: u32,
        /// Codomain dimension
        codomain_dim: u32,
        /// Element order
        order: u32,
        /// Whether the mesh is triangular
        triangular_mesh: bool,
    }

    unsafe extern "C++" {
        include!("pointsource_pde/include/dolfinx_access/nanobind_helpers.h");
        include!("pointsource_pde/include/dolfinx_access/function.h");
        include!("pointsource_pde/include/dolfinx_access/minmax_p2.h");

        /// C++ side Dolfinx `Function<double>`
        type Function_f64;
        /// C-side python object.
        type PyObject;

        /// Find the cell containing `x` (in 1D, 2D, or 3D, following Fenics weirdness
        /// of hard-coding vectors to 3D).
        pub unsafe fn cell_Function_f64(f: *const Function_f64, x: &[f64; 3]) -> i32;

        /// Evaluate `f` at `x`.
        pub unsafe fn eval_Function_f64(f: *const Function_f64, x: &[f64; 3]) -> f64;

        /// Evaluates `f` at `x` (in 1D, 2D, or 3D, following Fenics weirdness
        /// of hard-coding vectors to 3D), assuming `x` to belong to `cell`.
        pub unsafe fn eval_Function_f64_cell(
            f: *const Function_f64,
            x: &[f64; 3],
            cell: i32,
        ) -> f64;

        /// Evaluate the function at 6 coordinates known to belong to given cell
        pub unsafe fn eval_Function_f64_cell_6(
            f: *const Function_f64,
            all_coords: &[f64; 18], // 6*3
            cell: i32,
        ) -> [f64; 6];

        /// Get information about a dolfinx `Function_f64`.
        pub unsafe fn info_Function_f64(f: *const Function_f64) -> FunctionInfo;

        /// Returns the minimum value of the Dolfinx `Function`,
        /// if supported (P2 elements on triangular mesh)
        pub unsafe fn min_Function_f64_p2(f: *const Function_f64) -> Result<CoordValuePair>;

        /// Returns the maximum value of the Dolfinx `Function`,
        /// if supported (P2 elements on triangular mesh)
        pub unsafe fn max_Function_f64_p2(f: *const Function_f64) -> Result<CoordValuePair>;

        /// Returns the minimum or maximum value of the Dolfinx `Function`,
        /// if supported (P2 elements on triangular mesh)
        pub unsafe fn minmax_Function_f64_p2(
            f: *const Function_f64,
            max: bool,
        ) -> Result<CoordValuePair>;

        /// Checks if a Dolfinx `Function` is supported (P2 elements on triangular mesh)
        pub unsafe fn check_Function_f64(o: *const PyObject) -> bool;

        /// Casts a Python object into a Dolfinx `Function` wrapper.
        pub unsafe fn cast_Function_f64(o: *const PyObject) -> Result<*const Function_f64>;

        /// Casts a Python object into a mutable Dolfinx `Function` wrapper.
        pub unsafe fn cast_mut_Function_f64(o: *mut PyObject) -> Result<*mut Function_f64>;

        /// Gets the 3D coordinates of the nodes of a triangular cell
        pub unsafe fn cell_coords_Function_f64_triangle(
            f: *const Function_f64,
            cell: i32,
        ) -> [f64; 9];

        /// Returns the local weights
        pub unsafe fn data_Function_f64<'a>(f: *const Function_f64) -> &'a [f64];

        /// Returns the local mutable weights
        pub unsafe fn data_mut_Function_f64<'a>(f: *mut Function_f64) -> &'a mut [f64];

        /// Check the compatibility of function spaces of two functions.
        pub unsafe fn check_compat_Function_f64(
            f: *const Function_f64,
            g: *const Function_f64,
        ) -> Result<()>;

        /// Create a new similar Dolfinx `Function`
        pub unsafe fn similar_Function_f64(f: *const Function_f64) -> *mut Function_f64;

        /// Clone the Dolfinx `Function`
        pub unsafe fn clone_Function_f64(f: *const Function_f64) -> *mut Function_f64;

        /// Wraps a Dolfinx `Function` into a Python object
        pub unsafe fn wrap_Function_f64(f: *mut Function_f64) -> *mut PyObject;

        /// Drop helper
        pub unsafe fn drop_Function_f64(f: *mut Function_f64);

        /// Scatter helper. Not truly supported.
        pub unsafe fn scatter_fwd_Function_f64(f: *mut Function_f64);
    }

    extern "Rust" {
        /// Rust-side helper for minimising P2 elements. This bypasses poorly documented basix,
        /// and does things at a low level directly.
        /// `f` is the function, `simplex_coords` the coorindates of the triangle
        /// see [`cell_coords_Function_f64_triangle`], `cell` is the cell number,
        /// and `max` indicates whether to maximise or minimise.
        unsafe fn minmax_dolfinx_p2_cell(
            f: *const Function_f64,
            simplex_coords: &[f64; 9],
            cell: i32,
            max: bool,
        ) -> CoordValuePair;
    }
}

use ffi::CoordValuePair;

/// Convert an array of 2D oordinates obtained from Fenics, which hard-codes weird 3D coordinates.
#[inline]
fn coords_from_fenics(&[x1, x2, _, y1, y2, _, z1, z2, _]: &[f64; 9]) -> [Loc<2, f64>; 3] {
    [Loc([x1, x2]), Loc([y1, y2]), Loc([z1, z2])]
}
// fn coords_from_fenics<const N: usize>(coords: &[[f64; 3]; N]) -> [Loc<F, 2>; N] {
//     coords.map(|[x1, x2, _]| Loc([x1, x2]))
// }

/// Convert an array of 2D oordinates to Fenics, which hard-codes weird 3D coordinates.
#[inline]
fn coords_to_fenics([Loc([x1, x2]), Loc([y1, y2]), Loc([z1, z2])]: [Loc<2, f64>; 3]) -> [f64; 9] {
    [x1, x2, 0.0, y1, y2, 0.0, z1, z2, 0.0]
}
// fn coords_to_fenics<const N: usize>(coords: [Loc<F, 2>; N]) -> [[f64; 3]; N] {
//     coords.map(|Loc([x1, x2])| [x1, x2, 0.0])
// }

/// Given a dolphinx `Function` `f`, as well as 6 coordinates on a simplex in its mesh
/// (the corners and midpoints), construct a new internal P2 model on `cell`.
/// If `max` is true, the model is negated compared to `f`.
unsafe fn model_dolfinx_p2_cell(
    f: *const ffi::Function_f64,
    simplex_coords: &[f64; 9],
    cell: i32,
    max: bool,
) -> (P2LocalModel<f64, 2, 3>, Simplex<f64, 2, 3>) {
    // Form simplex
    let simplex = Simplex(coords_from_fenics(simplex_coords));
    // Form an array of all coordinates
    let midpoints = simplex.midpoints();
    let mut all_coords_uninit: [MaybeUninit<f64>; 6 * 3] = [const { MaybeUninit::uninit() }; 6 * 3];
    all_coords_uninit
        .iter_mut()
        .zip(
            simplex_coords
                .iter()
                .copied()
                .chain(coords_to_fenics(midpoints).into_iter()),
        )
        .for_each(|(a, b)| {
            a.write(b);
        });
    let all_coords = MaybeUninit::array_assume_init(all_coords_uninit);
    // Evaluate the function, and invert the signs of maximising
    let mut values = ffi::eval_Function_f64_cell_6(f, &all_coords, cell);
    if max {
        values.iter_mut().for_each(|v| *v = -*v);
    }
    // Form local model
    let [ref n0, ref n1, ref n2] = simplex.0;
    let [ref n01, ref n12, ref n20] = midpoints;
    let nodes = [*n0, *n1, *n2, *n01, *n12, *n20];
    let model = P2LocalModel::<f64, 2, 3>::new(&nodes, &values);

    return (model, simplex);
}

/// This is a helper routine for `minmax_dolfinx_p2`. It reconstructs a P2
/// model on the cell by evaluating f at 6 points (the simplex corners given
/// in `simplex_doords`, as well as edge midpoints). Then it minimises this
/// model. The result should be accurate, given that `f` is assumed to have a
/// P2 model on the simplex.
pub unsafe fn minmax_dolfinx_p2_cell(
    f: *const ffi::Function_f64,
    simplex_coords: &[f64; 9],
    cell: i32,
    max: bool,
) -> ffi::CoordValuePair {
    let (model, simplex) = model_dolfinx_p2_cell(f, simplex_coords, cell, max);
    let (x, v) = model.minimise(&simplex);
    CoordValuePair { x: x.into(), v }
}

pub const NAME: &str = "dolfinx_access";

unsafe extern "C" {
    pub fn PyInit__dolfinx_access() -> *mut pyo3::ffi::PyObject;
}

/// Register the module in Python
pub fn register_python_ffi() -> DynResult<()> {
    let cname = CString::new(NAME).unwrap();
    // We need to use into_raw() to transfer ownership; otherwise this stops working.
    if unsafe { PyImport_AppendInittab(cname.into_raw(), Some(PyInit__dolfinx_access)) } != 0 {
        bail!("Failed to add embedded nanobind dolfinx_access module");
    }
    Ok(())
}