gemlab 2.0.0

use super::get_interp_matrix;
use crate::integ::Gauss;
use crate::shapes::GeoClass;
use crate::shapes::Scratchpad;
use crate::StrError;
use russell_lab::{mat_inverse, mat_pseudo_inverse, Matrix};

/// Calculates the extrapolation matrix (integration points to nodes)
///
/// ```text
/// u_nodal =   E    u_points
///  (nv)    (nv,np)   (np)
/// ```
///
/// Where `nv` is the number of vertices (nodes) and `np` is the number of integration points (`ngauss`).
///
/// # Input
///
/// * `pad` -- The scratchpad
/// * `integ_points` -- Integration points' constants (ngauss)
///
/// # Output
///
/// * `E` -- The (nnode,ngauss) extrapolation matrix; aka "inverse" of the
///    interpolation matrix `P` calculated by [get_interp_matrix()]. See note below
///    regarding the inversion problem.
///
/// This function returns the "inverse" of the interpolation matrix `P`; however, the
/// inverse is only possible if `ngauss == nnode`. If there are more interpolation
/// points than nodes (`ngauss > node`), the problem is over-determined, and we use
/// the pseudo-inverse instead `P⁺`. If there are fewer interpolation points than nodes
/// (`ngauss < nnode`), the problem is under-determined and may even yield
/// spurious results. To avoid such results and minimize the errors, we adopt the method
/// proposed by Durand and Farias in Reference #1. In this case, a correction matrix `A` is
/// applied to the pseudo-inverse, and a translation `ξ B` is applied to the result.
/// An exception exists: with a single integration point, the pseudo-inverse is returned
/// without any correction. In summary:
///
/// ```text
///     ⎧ P⁻¹         if np = nv (exact)
/// E = ⎨ P⁺          if np > nv (over-determined)
///     │ P⁺          if np = 1 < nv (exception case)
///     ⎩ P⁺ A + ξ B  if np < nv (under-determined)
/// ```
///
/// where `ξ` is a matrix with the reference (natural) coordinates of all nodes (vertices),
/// augmented by a column of ones; see Eq. (32)-(34) of Ref #1. Let us define `hat_ξ`
/// as a matrix with the reference coordinates all integration points (sampling points),
/// augmented by a column of ones; see Eq. (32)-(34). Then we have that:
///
/// ```text
/// B = pseudo_inverse(hat_ξ)
/// A = I - hat_ξ B
/// ```
///
/// where `I` is an identity matrix with ngausss rows and columns.
///
/// With `nd` being the number of geometry dimensions plus 1, we have that:
///
/// ```text
/// P : np × nv     hat_ξ: np × nd
/// P⁺: nv × np     B    : nd × np
/// ξ : nv × nd     I    : np × np,   
/// ```
///
/// Thus
///
/// ```text
///    A   ←   I   - hat_ξ  *   B
/// (np,np) (np,np) (np,nd)  (nd,np)
///
///    E   ←   P⁺  *   A   +   ξ   *   B
/// (nv,np) (nv,np) (np,np) (nv,nd) (nd,np)
/// ```
///
/// # Reference
///
/// 1. Durand R and Farias MM (2014) A local extrapolation method for finite elements,
///    Advances in Engineering Software, 67:1-9 <https://doi.org/10.1016/j.advengsoft.2013.07.002>
///
/// # Examples
///
/// ```
/// use gemlab::integ::Gauss;
/// use gemlab::recovery::{get_extrap_matrix, get_interp_matrix};
/// use gemlab::shapes::{GeoKind, Scratchpad};
/// use gemlab::StrError;
/// use russell_lab::{mat_vec_mul, vec_approx_eq, Vector};
///
/// fn main() -> Result<(), StrError> {
///     //  6 2
///     //  5 | `.    * indicates the
///     //  4 | * `.    location of ips
///     //  3 |     `.
///     //  2 |       `.
///     //  1 | *     * `.
///     //  0 0-----------1
///     //    0 1 2 3 4 5 6
///
///     let space_ndim = 2;
///     let mut pad = Scratchpad::new(space_ndim, GeoKind::Tri3)?;
///     pad.set_xx(0, 0, 0.0);
///     pad.set_xx(0, 1, 0.0);
///     pad.set_xx(1, 0, 6.0);
///     pad.set_xx(1, 1, 0.0);
///     pad.set_xx(2, 0, 0.0);
///     pad.set_xx(2, 1, 6.0);
///
///     // nodal values
///     let u_nodal = Vector::from(&[1.0, 2.0, 3.0]);
///
///     // interpolated values
///     let gauss = Gauss::new_sized(pad.kind.class(), 3)?;
///     let mut u_points = Vector::new(gauss.npoint());
///     let pp = get_interp_matrix(&mut pad, &gauss);
///     mat_vec_mul(&mut u_points, 1.0, &pp, &u_nodal)?;
///
///     // extrapolated values (recovered)
///     let nnode = pad.xxt.dims().1;
///     let mut u_nodal_rec = Vector::new(nnode);
///     let ee = get_extrap_matrix(&mut pad, &gauss)?;
///     mat_vec_mul(&mut u_nodal_rec, 1.0, &ee, &u_points)?;
///
///     // check
///     vec_approx_eq(&u_nodal_rec, &u_nodal, 1e-14);
///     Ok(())
/// }
/// ```
pub fn get_extrap_matrix(pad: &mut Scratchpad, gauss: &Gauss) -> Result<Matrix, StrError> {
    // constants
    let (nv, geo_ndim) = pad.deriv.dims();
    let np = gauss.npoint();

    // calculate interpolation matrix P
    let mut pp = get_interp_matrix(pad, gauss);

    // allocate output
    let mut ee = Matrix::new(nv, np);

    // handle cases
    if np == nv {
        // E ← P⁻¹
        mat_inverse(&mut ee, &pp)?;
    } else if np > nv {
        // E ← P⁺
        mat_pseudo_inverse(&mut ee, &mut pp)?;
    } else if np == 1 {
        // E ← P⁺
        mat_pseudo_inverse(&mut ee, &mut pp)?;
    } else {
        // P⁺ ← pseudo-inverse of interpolation matrix
        let mut ppi = Matrix::new(nv, np);
        mat_pseudo_inverse(&mut ppi, &mut pp)?;

        // Bᵀ ← transpose(pinv(hat_ξ))
        let bbt = get_tr_pinv_hxi(pad.kind.class(), np);

        //    A   ←   I   - hat_ξ *   B
        // (np,np) (np,np) (np,nd) (nd,np)
        let mut aa = Matrix::new(np, np);
        for i in 0..np {
            for j in 0..np {
                let mut sum = 0.0;
                for k in 0..geo_ndim {
                    sum += gauss.coords(i)[k] * bbt[j][k];
                }
                sum += 1.0 * bbt[j][geo_ndim]; // 1.0 = gauss.coords(i)[geo_ndim]
                let d_ij = if i == j { 1.0 } else { 0.0 };
                aa.set(i, j, d_ij - sum);
            }
        }

        //    E   ←   ξ   *   B
        // (nv,np) (nv,nd) (nd,np)
        for i in 0..nv {
            let rci = pad.kind.reference_coords(i);
            for j in 0..np {
                for k in 0..geo_ndim {
                    ee.add(i, j, rci[k] * bbt[j][k]);
                }
                ee.add(i, j, 1.0 * bbt[j][geo_ndim]); // 1.0 = rci[geo_ndim]
            }
        }

        //    E   +=   P⁺  *   A
        // (nv,np)  (nv,np) (np,np)
        for i in 0..nv {
            for j in 0..np {
                for k in 0..np {
                    ee.add(i, j, ppi.get(i, k) * aa.get(k, j));
                }
            }
        }
    }
    return Ok(ee);
}

/// Returns the transpose of the pseudo-inverse of the reference coordinates matrix
///
/// # Output
///
/// ```text
/// TR_PINV_HXI ← transpose(PINV_HXI) ← transpose(pinv(hat_ξ))
/// ```
///
/// where `PINV_HXI ← pinv(hat_ξ)` and `HXI ← hat_ξ` is a matrix with the reference
/// (natural) coordinates of the integration points, augmented by a column of ones.
/// Here, `pinv` means the pseudo-inverse. See Reference #1.
///
/// The reason for returning the transpose is due to the way static arrays are
/// handled in Rust. With the transpose, we can share a references to the static
/// data from any GeoClass. To do so, we also save the transpose matrix with
/// 4 columns, even if geo_ndim < 3.
///
/// Note that the columns of the resulting matrix are padded with zeros if geo_ndim < 3.
///
/// # Reference
///
/// 1. Durand R and Farias MM (2014) A local extrapolation method for finite elements,
///    Advances in Engineering Software, 67:1-9 <https://doi.org/10.1016/j.advengsoft.2013.07.002>
fn get_tr_pinv_hxi(class: GeoClass, ngauss: usize) -> &'static [[f64; 4]] {
    match class {
        // Lin
        GeoClass::Lin => match ngauss {
            1 => &TR_PINV_HXI_LIN_LEGENDRE_1,
            2 => &TR_PINV_HXI_LIN_LEGENDRE_2,
            3 => &TR_PINV_HXI_LIN_LEGENDRE_3,
            4 => &TR_PINV_HXI_LIN_LEGENDRE_4,
            5 => &TR_PINV_HXI_LIN_LEGENDRE_5,
            _ => unreachable!("requested number of integration points is not available for Lin class"),
        },
        // Tri
        GeoClass::Tri => match ngauss {
            1 => &TR_PINV_HXI_TRI_INTERNAL_1,
            3 => &TR_PINV_HXI_TRI_INTERNAL_3,
            4 => &TR_PINV_HXI_TRI_INTERNAL_4,
            6 => &TR_PINV_HXI_TRI_FELIPPA_6,
            7 => &TR_PINV_HXI_TRI_FELIPPA_7,
            12 => &TR_PINV_HXI_TRI_INTERNAL_12,
            16 => &TR_PINV_HXI_TRI_INTERNAL_16,
            _ => unreachable!("requested number of integration points is not available for Tri class"),
        },
        // Qua
        GeoClass::Qua => match ngauss {
            1 => &TR_PINV_HXI_QUA_LEGENDRE_1,
            4 => &TR_PINV_HXI_QUA_LEGENDRE_4,
            9 => &TR_PINV_HXI_QUA_LEGENDRE_9,
            16 => &TR_PINV_HXI_QUA_LEGENDRE_16,
            _ => unreachable!("requested number of integration points is not available for Qua class"),
        },
        // Tet
        GeoClass::Tet => match ngauss {
            1 => &TR_PINV_HXI_TET_INTERNAL_1,
            4 => &TR_PINV_HXI_TET_INTERNAL_4,
            5 => &TR_PINV_HXI_TET_INTERNAL_5,
            8 => &TR_PINV_HXI_TET_FELIPPA_8,
            14 => &TR_PINV_HXI_TET_FELIPPA_14,
            15 => &TR_PINV_HXI_TET_FELIPPA_15,
            24 => &TR_PINV_HXI_TET_FELIPPA_24,
            _ => unreachable!("requested number of integration points is not available for Tet class"),
        },
        // Hex
        GeoClass::Hex => match ngauss {
            6 => &TR_PINV_HXI_HEX_IRONS_6,
            8 => &TR_PINV_HXI_HEX_LEGENDRE_8,
            14 => &TR_PINV_HXI_HEX_IRONS_14,
            27 => &TR_PINV_HXI_HEX_LEGENDRE_27,
            64 => &TR_PINV_HXI_HEX_LEGENDRE_64,
            _ => unreachable!("requested number of integration points is not available for Hex class"),
        },
    }
}

// -----------------------------------------------------------------------
// -- LIN ----------------------------------------------------------------
// -----------------------------------------------------------------------

// TR_PINV_HXI_LIN_LEGENDRE_1: nrow = 1, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_LIN_LEGENDRE_1: [[f64; 4]; 1] = [
    [  0.000000000000000E+00,  1.000000000000000E+00,  0.000000000000000E+00,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_LIN_LEGENDRE_2: nrow = 2, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_LIN_LEGENDRE_2: [[f64; 4]; 2] = [
    [ -8.660254037844387E-01,  5.000000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  8.660254037844387E-01,  5.000000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_LIN_LEGENDRE_3: nrow = 3, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_LIN_LEGENDRE_3: [[f64; 4]; 3] = [
    [ -6.454972243679027E-01,  3.333333333333334E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  0.000000000000000E+00,  3.333333333333334E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  6.454972243679027E-01,  3.333333333333334E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_LIN_LEGENDRE_4: nrow = 4, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_LIN_LEGENDRE_4: [[f64; 4]; 4] = [
    [ -5.023295150965307E-01,  2.500000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [ -1.983222754244993E-01,  2.499999999999999E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  1.983222754244995E-01,  2.500000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  5.023295150965307E-01,  2.500000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_LIN_LEGENDRE_5: nrow = 5, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_LIN_LEGENDRE_5: [[f64; 4]; 5] = [
    [ -4.077809306723987E-01,  1.999999999999999E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [ -2.423111895475574E-01,  2.000000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  0.000000000000000E+00,  2.000000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  2.423111895475574E-01,  2.000000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
    [  4.077809306723987E-01,  2.000000000000000E-01,  0.000000000000000E+00,  0.000000000000000E+00 ],
];

// -----------------------------------------------------------------------
// -- TRI ----------------------------------------------------------------
// -----------------------------------------------------------------------

// TR_PINV_HXI_TRI_INTERNAL_1: nrow = 1, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TRI_INTERNAL_1: [[f64; 4]; 1] = [
    [  2.727272727272727E-01,  2.727272727272728E-01,  8.181818181818182E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_TRI_INTERNAL_3: nrow = 3, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TRI_INTERNAL_3: [[f64; 4]; 3] = [
    [ -2.000000000000000E+00, -1.999999999999999E+00,  1.666666666666666E+00,  0.000000000000000E+00 ],
    [  2.000000000000000E+00, -9.992007221626409E-16, -3.333333333333329E-01,  0.000000000000000E+00 ],
    [ -5.551115123125783E-16,  2.000000000000000E+00, -3.333333333333331E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_TRI_INTERNAL_4: nrow = 4, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TRI_INTERNAL_4: [[f64; 4]; 4] = [
    [  3.053113317719180E-16, -4.440892098500626E-16,  2.500000000000000E-01,  0.000000000000000E+00 ],
    [ -2.500000000000000E+00, -2.499999999999999E+00,  1.916666666666667E+00,  0.000000000000000E+00 ],
    [  2.500000000000000E+00,  0.000000000000000E+00, -5.833333333333335E-01,  0.000000000000000E+00 ],
    [ -2.220446049250313E-16,  2.500000000000000E+00, -5.833333333333333E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_TRI_FELIPPA_6: nrow = 6, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TRI_FELIPPA_6: [[f64; 4]; 6] = [
    [  5.277531727779488E-01,  5.277531727779489E-01, -1.851687818519657E-01,  0.000000000000000E+00 ],
    [ -5.277531727779498E-01, -5.551115123125783E-16,  3.425843909259835E-01,  0.000000000000000E+00 ],
    [  1.110223024625157E-16, -5.277531727779492E-01,  3.425843909259830E-01,  0.000000000000000E+00 ],
    [ -1.132956608748925E+00, -1.132956608748925E+00,  9.219710724992830E-01,  0.000000000000000E+00 ],
    [  1.132956608748925E+00,  1.110223024625157E-16, -2.109855362496415E-01,  0.000000000000000E+00 ],
    [  3.330669073875470E-16,  1.132956608748925E+00, -2.109855362496416E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_TRI_FELIPPA_7: nrow = 7, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TRI_FELIPPA_7: [[f64; 4]; 7] = [
    [ -1.065965106982872E+00, -1.065965106982872E+00,  8.535005475123904E-01,  0.000000000000000E+00 ],
    [  1.065965106982872E+00,  1.665334536937735E-16, -2.124645594704811E-01,  0.000000000000000E+00 ],
    [ -1.054711873393899E-15,  1.065965106982872E+00, -2.124645594704807E-01,  0.000000000000000E+00 ],
    [  6.284651069828712E-01,  6.284651069828718E-01, -2.761195951314379E-01,  0.000000000000000E+00 ],
    [ -6.284651069828718E-01,  4.996003610813204E-16,  3.523455118514330E-01,  0.000000000000000E+00 ],
    [  8.881784197001252E-16, -6.284651069828719E-01,  3.523455118514331E-01,  0.000000000000000E+00 ],
    [ -6.938893903907228E-18, -6.938893903907228E-18,  1.428571428571429E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_TRI_INTERNAL_12: nrow = 12, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TRI_INTERNAL_12: [[f64; 4]; 12] = [
    [  6.571568458494592E-01,  5.551115123125783E-16, -1.357189486164864E-01,  0.000000000000000E+00 ],
    [  5.551115123125783E-17,  6.571568458494583E-01, -1.357189486164860E-01,  0.000000000000000E+00 ],
    [ -6.571568458494572E-01, -6.571568458494574E-01,  5.214378972329713E-01,  0.000000000000000E+00 ],
    [  2.043772502523998E-01, -2.359223927328458E-16,  1.520758324920021E-02,  0.000000000000000E+00 ],
    [ -2.081668171172169E-16,  2.043772502523998E-01,  1.520758324920021E-02,  0.000000000000000E+00 ],
    [ -2.043772502524000E-01, -2.043772502524001E-01,  2.195848335016000E-01,  0.000000000000000E+00 ],
    [ -4.728527911487811E-01, -2.643678601563865E-01,  3.290735504350557E-01,  0.000000000000000E+00 ],
    [ -2.643678601563865E-01, -4.728527911487813E-01,  3.290735504350558E-01,  0.000000000000000E+00 ],
    [ -2.084849309923949E-01,  2.643678601563862E-01,  6.470569027866957E-02,  0.000000000000000E+00 ],
    [  2.084849309923942E-01,  4.728527911487809E-01, -1.437792407137249E-01,  0.000000000000000E+00 ],
    [  2.643678601563864E-01, -2.084849309923949E-01,  6.470569027866956E-02,  0.000000000000000E+00 ],
    [  4.728527911487808E-01,  2.084849309923944E-01, -1.437792407137249E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_TRI_INTERNAL_16: nrow = 16, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TRI_INTERNAL_16: [[f64; 4]; 16] = [
    [ -7.702172233337023E-16, -7.632783294297951E-16,  6.250000000000060E-02,  0.000000000000000E+00 ],
    [ -1.987854631147938E-01,  8.326672684688674E-16,  1.287618210382643E-01,  0.000000000000000E+00 ],
    [  8.326672684688674E-17, -1.987854631147941E-01,  1.287618210382647E-01,  0.000000000000000E+00 ],
    [  1.987854631147946E-01,  1.987854631147943E-01, -7.002364207652959E-02,  0.000000000000000E+00 ],
    [  2.568697489842479E-01, -9.714451465470120E-17, -2.312324966141592E-02,  0.000000000000000E+00 ],
    [  1.249000902703301E-16,  2.568697489842479E-01, -2.312324966141593E-02,  0.000000000000000E+00 ],
    [ -2.568697489842480E-01, -2.568697489842476E-01,  2.337464993228319E-01,  0.000000000000000E+00 ],
    [  4.462853235057508E-01, -8.326672684688674E-17, -8.626177450191685E-02,  0.000000000000000E+00 ],
    [  1.110223024625157E-16,  4.462853235057505E-01, -8.626177450191681E-02,  0.000000000000000E+00 ],
    [ -4.462853235057508E-01, -4.462853235057501E-01,  3.600235490038337E-01,  0.000000000000000E+00 ],
    [ -3.788128101165515E-01, -2.448164475851215E-01,  2.703764192338910E-01,  0.000000000000000E+00 ],
    [  2.448164475851219E-01, -1.339963625314294E-01,  2.555997164876920E-02,  0.000000000000000E+00 ],
    [  1.339963625314296E-01,  3.788128101165512E-01, -1.084363908826602E-01,  0.000000000000000E+00 ],
    [ -1.339963625314294E-01,  2.448164475851219E-01,  2.555997164876921E-02,  0.000000000000000E+00 ],
    [  3.788128101165515E-01,  1.339963625314294E-01, -1.084363908826603E-01,  0.000000000000000E+00 ],
    [ -2.448164475851220E-01, -3.788128101165510E-01,  2.703764192338910E-01,  0.000000000000000E+00 ],
];

// -----------------------------------------------------------------------
// -- QUA ----------------------------------------------------------------
// -----------------------------------------------------------------------

// TR_PINV_HXI_QUA_LEGENDRE_1: nrow = 1, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_QUA_LEGENDRE_1: [[f64; 4]; 1] = [
    [  0.000000000000000E+00,  0.000000000000000E+00,  1.000000000000000E+00,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_QUA_LEGENDRE_4: nrow = 4, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_QUA_LEGENDRE_4: [[f64; 4]; 4] = [
    [ -4.330127018922194E-01, -4.330127018922194E-01,  2.500000000000001E-01,  0.000000000000000E+00 ],
    [  4.330127018922194E-01, -4.330127018922194E-01,  2.500000000000001E-01,  0.000000000000000E+00 ],
    [ -4.330127018922194E-01,  4.330127018922192E-01,  2.500000000000001E-01,  0.000000000000000E+00 ],
    [  4.330127018922194E-01,  4.330127018922192E-01,  2.500000000000001E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_QUA_LEGENDRE_9: nrow = 9, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_QUA_LEGENDRE_9: [[f64; 4]; 9] = [
    [ -2.151657414559677E-01, -2.151657414559677E-01,  1.111111111111112E-01,  0.000000000000000E+00 ],
    [  0.000000000000000E+00, -2.151657414559677E-01,  1.111111111111111E-01,  0.000000000000000E+00 ],
    [  2.151657414559676E-01, -2.151657414559676E-01,  1.111111111111112E-01,  0.000000000000000E+00 ],
    [ -2.151657414559676E-01, -2.159566345137073E-17,  1.111111111111111E-01,  0.000000000000000E+00 ],
    [  0.000000000000000E+00,  0.000000000000000E+00,  1.111111111111111E-01,  0.000000000000000E+00 ],
    [  2.151657414559676E-01,  2.159566345137073E-17,  1.111111111111111E-01,  0.000000000000000E+00 ],
    [ -2.151657414559676E-01,  2.151657414559676E-01,  1.111111111111111E-01,  0.000000000000000E+00 ],
    [  0.000000000000000E+00,  2.151657414559676E-01,  1.111111111111111E-01,  0.000000000000000E+00 ],
    [  2.151657414559676E-01,  2.151657414559677E-01,  1.111111111111111E-01,  0.000000000000000E+00 ],
];

// TR_PINV_HXI_QUA_LEGENDRE_16: nrow = 16, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_QUA_LEGENDRE_16: [[f64; 4]; 16] = [
    [ -1.255823787741326E-01, -1.255823787741326E-01,  6.249999999999997E-02,  0.000000000000000E+00 ],
    [ -4.958056885612486E-02, -1.255823787741327E-01,  6.249999999999997E-02,  0.000000000000000E+00 ],
    [  4.958056885612486E-02, -1.255823787741326E-01,  6.249999999999996E-02,  0.000000000000000E+00 ],
    [  1.255823787741326E-01, -1.255823787741327E-01,  6.249999999999998E-02,  0.000000000000000E+00 ],
    [ -1.255823787741326E-01, -4.958056885612488E-02,  6.249999999999999E-02,  0.000000000000000E+00 ],
    [ -4.958056885612486E-02, -4.958056885612488E-02,  6.249999999999998E-02,  0.000000000000000E+00 ],
    [  4.958056885612486E-02, -4.958056885612488E-02,  6.249999999999998E-02,  0.000000000000000E+00 ],
    [  1.255823787741326E-01, -4.958056885612488E-02,  6.249999999999999E-02,  0.000000000000000E+00 ],
    [ -1.255823787741326E-01,  4.958056885612484E-02,  6.249999999999997E-02,  0.000000000000000E+00 ],
    [ -4.958056885612486E-02,  4.958056885612485E-02,  6.249999999999997E-02,  0.000000000000000E+00 ],
    [  4.958056885612486E-02,  4.958056885612485E-02,  6.249999999999995E-02,  0.000000000000000E+00 ],
    [  1.255823787741326E-01,  4.958056885612484E-02,  6.249999999999997E-02,  0.000000000000000E+00 ],
    [ -1.255823787741326E-01,  1.255823787741326E-01,  6.249999999999997E-02,  0.000000000000000E+00 ],
    [ -4.958056885612486E-02,  1.255823787741326E-01,  6.249999999999998E-02,  0.000000000000000E+00 ],
    [  4.958056885612486E-02,  1.255823787741326E-01,  6.249999999999996E-02,  0.000000000000000E+00 ],
    [  1.255823787741326E-01,  1.255823787741326E-01,  6.249999999999997E-02,  0.000000000000000E+00 ],
];

// -----------------------------------------------------------------------
// -- TET ----------------------------------------------------------------
// -----------------------------------------------------------------------

// TR_PINV_HXI_TET_INTERNAL_1: nrow = 1, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TET_INTERNAL_1: [[f64; 4]; 1] = [
    [  2.105263157894737E-01,  2.105263157894737E-01,  2.105263157894737E-01,  8.421052631578949E-01 ],
];

// TR_PINV_HXI_TET_INTERNAL_4: nrow = 4, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TET_INTERNAL_4: [[f64; 4]; 4] = [
    [ -2.236067977499792E+00, -2.236067977499792E+00, -2.236067977499791E+00,  1.927050983124843E+00 ],
    [  2.236067977499790E+00, -1.110223024625157E-16,  4.440892098500626E-16, -3.090169943749476E-01 ],
    [  1.332267629550188E-15,  2.236067977499791E+00,  9.992007221626409E-16, -3.090169943749480E-01 ],
    [  6.661338147750939E-16,  6.661338147750939E-16,  2.236067977499790E+00, -3.090169943749475E-01 ],
];

// TR_PINV_HXI_TET_INTERNAL_5: nrow = 5, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TET_INTERNAL_5: [[f64; 4]; 5] = [
    [  5.134781488891349E-16,  4.996003610813204E-16,  5.412337245047638E-16,  1.999999999999996E-01 ],
    [ -3.000000000000000E+00, -3.000000000000000E+00, -3.000000000000000E+00,  2.450000000000000E+00 ],
    [  1.110223024625157E-16, -2.220446049250313E-16,  3.000000000000000E+00, -5.499999999999998E-01 ],
    [ -2.220446049250313E-16,  3.000000000000000E+00, -6.661338147750939E-16, -5.499999999999996E-01 ],
    [  3.000000000000000E+00, -3.330669073875470E-16, -4.440892098500626E-16, -5.499999999999997E-01 ],
];

// TR_PINV_HXI_TET_FELIPPA_8: nrow = 8, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TET_FELIPPA_8: [[f64; 4]; 8] = [
    [  7.348347322268437E-01,  7.348347322268440E-01,  7.348347322268439E-01, -4.261260491701330E-01 ],
    [ -7.348347322268431E-01,  2.775557561562891E-16,  3.885780586188048E-16,  3.087086830567106E-01 ],
    [ -2.775557561562891E-16, -7.348347322268443E-01, -3.330669073875470E-16,  3.087086830567113E-01 ],
    [  0.000000000000000E+00, -1.665334536937735E-16, -7.348347322268436E-01,  3.087086830567109E-01 ],
    [ -1.346702273007044E+00, -1.346702273007045E+00, -1.346702273007044E+00,  1.135026704755283E+00 ],
    [  1.346702273007045E+00,  3.330669073875470E-16, -1.110223024625157E-16, -2.116755682517612E-01 ],
    [ -4.996003610813204E-16,  1.346702273007044E+00, -3.330669073875470E-16, -2.116755682517609E-01 ],
    [ -3.885780586188048E-16, -8.881784197001252E-16,  1.346702273007044E+00, -2.116755682517606E-01 ],
];

// TR_PINV_HXI_TET_FELIPPA_14: nrow = 14, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TET_FELIPPA_14: [[f64; 4]; 14] = [
    [ -7.967021616808494E-01, -7.967021616808487E-01, -7.967021616808488E-01,  6.689551926892084E-01 ],
    [  7.967021616808486E-01, -4.718447854656915E-16, -5.273559366969494E-16, -1.277469689916405E-01 ],
    [  7.771561172376096E-16,  7.967021616808496E-01,  5.828670879282072E-16, -1.277469689916413E-01 ],
    [ -2.775557561562891E-16, -7.771561172376096E-16,  7.967021616808483E-01, -1.277469689916405E-01 ],
    [  3.084476565084738E-01,  3.084476565084736E-01,  3.084476565084736E-01, -1.599071709527839E-01 ],
    [ -3.084476565084734E-01, -2.220446049250313E-16, -1.526556658859590E-16,  1.485404855556900E-01 ],
    [ -5.689893001203927E-16, -3.084476565084737E-01, -4.440892098500626E-16,  1.485404855556902E-01 ],
    [  1.526556658859590E-16,  3.191891195797325E-16, -3.084476565084729E-01,  1.485404855556896E-01 ],
    [  7.771561172376096E-16,  5.179884281153202E-01,  5.179884281153203E-01, -1.875656426290890E-01 ],
    [  5.179884281153202E-01, -5.551115123125783E-17,  5.179884281153199E-01, -1.875656426290887E-01 ],
    [  5.179884281153215E-01,  5.179884281153214E-01,  1.443289932012704E-15, -1.875656426290898E-01 ],
    [ -5.179884281153214E-01, -5.179884281153213E-01, -1.332267629550188E-15,  3.304227854862326E-01 ],
    [ -5.179884281153200E-01,  2.775557561562891E-16, -5.179884281153196E-01,  3.304227854862315E-01 ],
    [ -5.551115123125783E-16, -5.179884281153200E-01, -5.179884281153200E-01,  3.304227854862317E-01 ],
];

// TR_PINV_HXI_TET_FELIPPA_15: nrow = 15, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TET_FELIPPA_15: [[f64; 4]; 15] = [
    [ -8.130014854074423E-01, -8.130014854074422E-01, -8.130014854074420E-01,  6.764177807222482E-01 ],
    [  8.130014854074428E-01,  4.996003610813204E-16,  2.775557561562891E-16, -1.365837046851943E-01 ],
    [ -6.661338147750939E-16,  8.130014854074419E-01, -3.885780586188048E-16, -1.365837046851935E-01 ],
    [ -3.330669073875470E-16, -2.775557561562891E-16,  8.130014854074420E-01, -1.365837046851937E-01 ],
    [  3.590629006277350E-01,  3.590629006277352E-01,  3.590629006277353E-01, -2.026305088041347E-01 ],
    [ -3.590629006277354E-01,  4.163336342344337E-16,  3.608224830031759E-16,  1.564323918236003E-01 ],
    [  1.942890293094024E-16, -3.590629006277357E-01,  1.942890293094024E-16,  1.564323918236004E-01 ],
    [  2.775557561562891E-17,  1.110223024625157E-16, -3.590629006277356E-01,  1.564323918236005E-01 ],
    [  2.220446049250313E-16, -4.981273640649506E-01, -4.981273640649505E-01,  3.157303486991418E-01 ],
    [ -4.981273640649502E-01,  6.106226635438361E-16, -4.981273640649501E-01,  3.157303486991416E-01 ],
    [ -4.981273640649498E-01, -4.981273640649499E-01,  7.216449660063518E-16,  3.157303486991413E-01 ],
    [  4.981273640649497E-01,  4.981273640649500E-01, -7.216449660063518E-16, -1.823970153658079E-01 ],
    [  4.981273640649501E-01, -7.216449660063518E-16,  4.981273640649502E-01, -1.823970153658082E-01 ],
    [ -4.440892098500626E-16,  4.981273640649506E-01,  4.981273640649506E-01, -1.823970153658085E-01 ],
    [  3.122502256758253E-17,  1.092875789865388E-16,  1.387778780781446E-16,  6.666666666666658E-02 ],
];

// TR_PINV_HXI_TET_FELIPPA_24: nrow = 24, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_TET_FELIPPA_24: [[f64; 4]; 24] = [
    [ -8.946419742957495E-02, -8.946419742957513E-02, -8.946419742957508E-02,  1.087648147388479E-01 ],
    [  8.946419742957515E-02,  9.020562075079397E-17,  6.591949208711867E-17,  1.930061730927279E-02 ],
    [ -6.938893903907228E-17,  8.946419742957507E-02, -1.457167719820518E-16,  1.930061730927288E-02 ],
    [  5.204170427930421E-17,  1.110223024625157E-16,  8.946419742957519E-02,  1.930061730927277E-02 ],
    [ -5.290594736629001E-01, -5.290594736629005E-01, -5.290594736629003E-01,  4.384612719138419E-01 ],
    [  5.290594736629005E-01,  2.775557561562891E-17,  2.775557561562891E-17, -9.059820174905855E-02 ],
    [ -9.159339953157541E-16,  5.290594736628997E-01, -1.026956297778270E-15, -9.059820174905792E-02 ],
    [  6.383782391594650E-16,  8.049116928532385E-16,  5.290594736629013E-01, -9.059820174905911E-02 ],
    [  1.828298372082145E-01,  1.828298372082146E-01,  1.828298372082145E-01, -9.545571123949428E-02 ],
    [ -1.828298372082147E-01, -6.938893903907228E-17, -2.775557561562891E-17,  8.737412596872032E-02 ],
    [  2.220446049250313E-16, -1.828298372082145E-01,  2.081668171172169E-16,  8.737412596872016E-02 ],
    [ -3.469446951953614E-16, -3.885780586188048E-16, -1.828298372082151E-01,  8.737412596872057E-02 ],
    [ -2.106201844895547E-01, -3.407906172208732E-01, -3.407906172208731E-01,  2.647170213994920E-01 ],
    [ -3.407906172208732E-01, -2.106201844895551E-01, -3.407906172208733E-01,  2.647170213994921E-01 ],
    [ -3.407906172208729E-01, -3.407906172208731E-01, -2.106201844895546E-01,  2.647170213994919E-01 ],
    [  3.407906172208731E-01,  1.301704327313181E-01, -2.220446049250313E-16, -7.607359582138114E-02 ],
    [  3.407906172208734E-01,  1.387778780781446E-16,  1.301704327313185E-01, -7.607359582138142E-02 ],
    [ -4.163336342344337E-16,  3.407906172208730E-01,  1.301704327313179E-01, -7.607359582138104E-02 ],
    [  2.106201844895549E-01, -1.301704327313183E-01, -1.301704327313183E-01,  5.409683690993707E-02 ],
    [ -1.301704327313188E-01,  2.106201844895546E-01, -1.301704327313189E-01,  5.409683690993741E-02 ],
    [ -1.301704327313179E-01, -1.301704327313179E-01,  2.106201844895555E-01,  5.409683690993669E-02 ],
    [  1.301704327313178E-01,  3.407906172208728E-01, -6.106226635438361E-16, -7.607359582138090E-02 ],
    [  1.301704327313188E-01,  4.996003610813204E-16,  3.407906172208738E-01, -7.607359582138166E-02 ],
    [  1.942890293094024E-16,  1.301704327313186E-01,  3.407906172208736E-01, -7.607359582138153E-02 ],
];

// -----------------------------------------------------------------------
// -- HEX ----------------------------------------------------------------
// -----------------------------------------------------------------------

// TR_PINV_HXI_HEX_IRONS_6: nrow = 6, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_HEX_IRONS_6: [[f64; 4]; 6] = [
    [ -4.999999999999998E-01,  0.000000000000000E+00,  0.000000000000000E+00,  1.666666666666666E-01 ],
    [  4.999999999999999E-01,  0.000000000000000E+00,  0.000000000000000E+00,  1.666666666666666E-01 ],
    [  0.000000000000000E+00, -4.999999999999999E-01,  0.000000000000000E+00,  1.666666666666667E-01 ],
    [  0.000000000000000E+00,  4.999999999999999E-01,  0.000000000000000E+00,  1.666666666666666E-01 ],
    [  0.000000000000000E+00,  0.000000000000000E+00, -4.999999999999999E-01,  1.666666666666668E-01 ],
    [  0.000000000000000E+00,  0.000000000000000E+00,  4.999999999999999E-01,  1.666666666666667E-01 ],
];

// TR_PINV_HXI_HEX_LEGENDRE_8: nrow = 8, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_HEX_LEGENDRE_8: [[f64; 4]; 8] = [
    [ -2.165063509461095E-01, -2.165063509461095E-01, -2.165063509461096E-01,  1.250000000000000E-01 ],
    [  2.165063509461096E-01, -2.165063509461097E-01, -2.165063509461097E-01,  1.250000000000000E-01 ],
    [ -2.165063509461096E-01,  2.165063509461095E-01, -2.165063509461096E-01,  1.250000000000000E-01 ],
    [  2.165063509461096E-01,  2.165063509461096E-01, -2.165063509461096E-01,  1.250000000000001E-01 ],
    [ -2.165063509461096E-01, -2.165063509461095E-01,  2.165063509461096E-01,  1.250000000000000E-01 ],
    [  2.165063509461096E-01, -2.165063509461094E-01,  2.165063509461095E-01,  1.250000000000000E-01 ],
    [ -2.165063509461096E-01,  2.165063509461097E-01,  2.165063509461097E-01,  1.250000000000000E-01 ],
    [  2.165063509461096E-01,  2.165063509461097E-01,  2.165063509461096E-01,  1.250000000000000E-01 ],
];

// TR_PINV_HXI_HEX_IRONS_14: nrow = 14, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_HEX_IRONS_14: [[f64; 4]; 14] = [
    [  1.355115585649603E-01,  2.290658341603791E-17,  0.000000000000000E+00,  7.142857142857142E-02 ],
    [ -1.355115585649603E-01,  2.290658341603791E-17,  0.000000000000000E+00,  7.142857142857142E-02 ],
    [  0.000000000000000E+00,  1.355115585649603E-01,  0.000000000000000E+00,  7.142857142857142E-02 ],
    [  0.000000000000000E+00, -1.355115585649603E-01,  0.000000000000000E+00,  7.142857142857144E-02 ],
    [  0.000000000000000E+00,  2.290658341603791E-17,  1.355115585649603E-01,  7.142857142857142E-02 ],
    [  0.000000000000000E+00,  2.290658341603791E-17, -1.355115585649603E-01,  7.142857142857142E-02 ],
    [  1.292052015020528E-01,  1.292052015020528E-01,  1.292052015020528E-01,  7.142857142857142E-02 ],
    [ -1.292052015020528E-01,  1.292052015020528E-01,  1.292052015020528E-01,  7.142857142857142E-02 ],
    [  1.292052015020528E-01, -1.292052015020528E-01,  1.292052015020528E-01,  7.142857142857145E-02 ],
    [ -1.292052015020528E-01, -1.292052015020527E-01,  1.292052015020528E-01,  7.142857142857142E-02 ],
    [  1.292052015020528E-01,  1.292052015020528E-01, -1.292052015020528E-01,  7.142857142857142E-02 ],
    [ -1.292052015020528E-01,  1.292052015020528E-01, -1.292052015020528E-01,  7.142857142857142E-02 ],
    [  1.292052015020528E-01, -1.292052015020527E-01, -1.292052015020528E-01,  7.142857142857142E-02 ],
    [ -1.292052015020528E-01, -1.292052015020527E-01, -1.292052015020528E-01,  7.142857142857142E-02 ],
];

// TR_PINV_HXI_HEX_LEGENDRE_27: nrow = 27, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_HEX_LEGENDRE_27: [[f64; 4]; 27] = [
    [ -7.172191381865585E-02, -7.172191381865588E-02, -7.172191381865589E-02,  3.703703703703706E-02 ],
    [  0.000000000000000E+00, -7.172191381865585E-02, -7.172191381865592E-02,  3.703703703703706E-02 ],
    [  7.172191381865589E-02, -7.172191381865591E-02, -7.172191381865592E-02,  3.703703703703706E-02 ],
    [ -7.172191381865589E-02,  2.081668171172169E-17, -7.172191381865589E-02,  3.703703703703708E-02 ],
    [  0.000000000000000E+00,  0.000000000000000E+00, -7.172191381865589E-02,  3.703703703703708E-02 ],
    [  7.172191381865589E-02,  1.387778780781446E-17, -7.172191381865589E-02,  3.703703703703707E-02 ],
    [ -7.172191381865589E-02,  7.172191381865591E-02, -7.172191381865589E-02,  3.703703703703709E-02 ],
    [  0.000000000000000E+00,  7.172191381865591E-02, -7.172191381865591E-02,  3.703703703703708E-02 ],
    [  7.172191381865589E-02,  7.172191381865591E-02, -7.172191381865589E-02,  3.703703703703709E-02 ],
    [ -7.172191381865589E-02, -7.172191381865589E-02, -2.775557561562891E-17,  3.703703703703703E-02 ],
    [  0.000000000000000E+00, -7.172191381865586E-02, -6.938893903907228E-18,  3.703703703703703E-02 ],
    [  7.172191381865589E-02, -7.172191381865591E-02, -1.387778780781446E-17,  3.703703703703703E-02 ],
    [ -7.172191381865589E-02,  2.659175150280490E-18,  1.082140139430786E-18,  3.703703703703706E-02 ],
    [  0.000000000000000E+00, -3.204251063349822E-19, -1.622456720000829E-18,  3.703703703703704E-02 ],
    [  7.172191381865589E-02,  2.593636424395281E-18,  1.291604810729439E-18,  3.703703703703705E-02 ],
    [ -7.172191381865589E-02,  7.172191381865589E-02,  2.081668171172169E-17,  3.703703703703706E-02 ],
    [  0.000000000000000E+00,  7.172191381865589E-02,  2.081668171172169E-17,  3.703703703703706E-02 ],
    [  7.172191381865589E-02,  7.172191381865589E-02,  1.387778780781446E-17,  3.703703703703706E-02 ],
    [ -7.172191381865589E-02, -7.172191381865591E-02,  7.172191381865589E-02,  3.703703703703701E-02 ],
    [  0.000000000000000E+00, -7.172191381865591E-02,  7.172191381865589E-02,  3.703703703703701E-02 ],
    [  7.172191381865589E-02, -7.172191381865591E-02,  7.172191381865589E-02,  3.703703703703701E-02 ],
    [ -7.172191381865589E-02,  0.000000000000000E+00,  7.172191381865589E-02,  3.703703703703703E-02 ],
    [  0.000000000000000E+00, -6.938893903907228E-18,  7.172191381865589E-02,  3.703703703703701E-02 ],
    [  7.172191381865589E-02, -6.938893903907228E-18,  7.172191381865592E-02,  3.703703703703703E-02 ],
    [ -7.172191381865589E-02,  7.172191381865588E-02,  7.172191381865589E-02,  3.703703703703703E-02 ],
    [  0.000000000000000E+00,  7.172191381865588E-02,  7.172191381865589E-02,  3.703703703703705E-02 ],
    [  7.172191381865589E-02,  7.172191381865588E-02,  7.172191381865589E-02,  3.703703703703703E-02 ],
];

// TR_PINV_HXI_HEX_LEGENDRE_64: nrow = 64, ncol = 4
#[rustfmt::skip]
const TR_PINV_HXI_HEX_LEGENDRE_64: [[f64; 4]; 64] = [
    [ -3.139559469353317E-02, -3.139559469353318E-02, -3.139559469353316E-02,  1.562500000000000E-02 ],
    [ -1.239514221403122E-02, -3.139559469353315E-02, -3.139559469353317E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02, -3.139559469353317E-02, -3.139559469353320E-02,  1.562500000000000E-02 ],
    [  3.139559469353316E-02, -3.139559469353317E-02, -3.139559469353317E-02,  1.562500000000001E-02 ],
    [ -3.139559469353316E-02, -1.239514221403122E-02, -3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02, -1.239514221403122E-02, -3.139559469353317E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02, -1.239514221403121E-02, -3.139559469353317E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02, -1.239514221403122E-02, -3.139559469353318E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02,  1.239514221403122E-02, -3.139559469353317E-02,  1.562500000000000E-02 ],
    [ -1.239514221403122E-02,  1.239514221403121E-02, -3.139559469353318E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02,  1.239514221403122E-02, -3.139559469353317E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02,  1.239514221403122E-02, -3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02,  3.139559469353317E-02, -3.139559469353317E-02,  1.562500000000000E-02 ],
    [ -1.239514221403122E-02,  3.139559469353317E-02, -3.139559469353317E-02,  1.562500000000000E-02 ],
    [  1.239514221403122E-02,  3.139559469353317E-02, -3.139559469353319E-02,  1.562500000000000E-02 ],
    [  3.139559469353316E-02,  3.139559469353317E-02, -3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02, -3.139559469353317E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02, -3.139559469353317E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02, -3.139559469353317E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02, -3.139559469353317E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02, -1.239514221403122E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02, -1.239514221403122E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02, -1.239514221403121E-02, -1.239514221403123E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02, -1.239514221403122E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02,  1.239514221403122E-02, -1.239514221403122E-02,  1.562500000000000E-02 ],
    [ -1.239514221403122E-02,  1.239514221403121E-02, -1.239514221403122E-02,  1.562500000000000E-02 ],
    [  1.239514221403122E-02,  1.239514221403122E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02,  1.239514221403122E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02,  3.139559469353317E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02,  3.139559469353317E-02, -1.239514221403122E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02,  3.139559469353317E-02, -1.239514221403122E-02,  1.562500000000000E-02 ],
    [  3.139559469353316E-02,  3.139559469353317E-02, -1.239514221403122E-02,  1.562500000000000E-02 ],
    [ -3.139559469353316E-02, -3.139559469353317E-02,  1.239514221403123E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02, -3.139559469353317E-02,  1.239514221403122E-02,  1.562500000000000E-02 ],
    [  1.239514221403122E-02, -3.139559469353317E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02, -3.139559469353317E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02, -1.239514221403122E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02, -1.239514221403122E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02, -1.239514221403121E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02, -1.239514221403122E-02,  1.239514221403122E-02,  1.562500000000000E-02 ],
    [ -3.139559469353316E-02,  1.239514221403122E-02,  1.239514221403122E-02,  1.562500000000000E-02 ],
    [ -1.239514221403122E-02,  1.239514221403121E-02,  1.239514221403123E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02,  1.239514221403122E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02,  1.239514221403122E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02,  3.139559469353317E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02,  3.139559469353317E-02,  1.239514221403123E-02,  1.562500000000000E-02 ],
    [  1.239514221403122E-02,  3.139559469353317E-02,  1.239514221403122E-02,  1.562500000000000E-02 ],
    [  3.139559469353316E-02,  3.139559469353317E-02,  1.239514221403122E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02, -3.139559469353317E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02, -3.139559469353317E-02,  3.139559469353319E-02,  1.562500000000000E-02 ],
    [  1.239514221403122E-02, -3.139559469353317E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02, -3.139559469353317E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02, -1.239514221403122E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -1.239514221403122E-02, -1.239514221403122E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02, -1.239514221403121E-02,  3.139559469353318E-02,  1.562500000000000E-02 ],
    [  3.139559469353316E-02, -1.239514221403122E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02,  1.239514221403122E-02,  3.139559469353318E-02,  1.562500000000000E-02 ],
    [ -1.239514221403122E-02,  1.239514221403121E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02,  1.239514221403122E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [  3.139559469353316E-02,  1.239514221403122E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
    [ -3.139559469353316E-02,  3.139559469353317E-02,  3.139559469353317E-02,  1.562500000000000E-02 ],
    [ -1.239514221403122E-02,  3.139559469353317E-02,  3.139559469353318E-02,  1.562499999999999E-02 ],
    [  1.239514221403122E-02,  3.139559469353317E-02,  3.139559469353317E-02,  1.562500000000000E-02 ],
    [  3.139559469353316E-02,  3.139559469353317E-02,  3.139559469353317E-02,  1.562499999999999E-02 ],
];

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use super::get_extrap_matrix;
    use crate::integ::Gauss;
    use crate::recovery::{get_interp_matrix, get_points_coords};
    use crate::shapes::{GeoClass, GeoKind, Scratchpad};
    use plotpy::{generate3d, Plot, Surface};
    use russell_lab::math::PI;
    use russell_lab::{
        approx_eq, mat_approx_eq, mat_pseudo_inverse, mat_vec_mul, vec_approx_eq, vec_inner, Matrix, Vector,
    };

    const SAVE_FIGURE: bool = false;

    fn gen_qua4(w: f64, h: f64, skew_angle: f64, save_fig: bool) -> Scratchpad {
        //  3-------------2         ξ₀   ξ₁
        //  |      ξ₁     |  node    r    s
        //  |      |      |     0 -1.0 -1.0
        //  |      +--ξ₀  |     1  1.0 -1.0
        //  |             |     2  1.0  1.0
        //  |             |     3 -1.0  1.0
        //  0-------------1
        let space_ndim = 2;
        let mut pad = Scratchpad::new(space_ndim, GeoKind::Qua4).unwrap();
        let dx = h * f64::tan(skew_angle);
        pad.set_xx(0, 0, 0.0);
        pad.set_xx(0, 1, 0.0);
        pad.set_xx(1, 0, w);
        pad.set_xx(1, 1, 0.0);
        pad.set_xx(2, 0, w + dx);
        pad.set_xx(2, 1, h);
        pad.set_xx(3, 0, 0.0 + dx);
        pad.set_xx(3, 1, h);
        assert!(pad.ok_xxt);
        if save_fig {
            pad.draw_shape_simple("/tmp/gemlab/gen_qua4.svg").unwrap();
        }
        pad
    }

    fn gen_qua8(w: f64, h: f64, skew_angle: f64, imprecision: f64, save_fig: bool) -> Scratchpad {
        //  3------6------2         ξ₀   ξ₁
        //  |      ξ₁     |  node    r    s
        //  |      |      |     0 -1.0 -1.0   4  0.0 -1.0
        //  7      *--ξ₀  5     1  1.0 -1.0   5  1.0  0.0
        //  |             |     2  1.0  1.0   6  0.0  1.0
        //  |             |     3 -1.0  1.0   7 -1.0  0.0
        //  0------4------1
        let space_ndim = 2;
        let mut pad = Scratchpad::new(space_ndim, GeoKind::Qua8).unwrap();
        let dx = h * f64::tan(skew_angle);
        pad.set_xx(0, 0, 0.0);
        pad.set_xx(0, 1, 0.0);
        pad.set_xx(1, 0, w);
        pad.set_xx(1, 1, 0.0);
        pad.set_xx(2, 0, w + dx);
        pad.set_xx(2, 1, h);
        pad.set_xx(3, 0, 0.0 + dx);
        pad.set_xx(3, 1, h);
        pad.set_xx(4, 0, w / 2.0);
        pad.set_xx(4, 1, 0.0 + imprecision);
        pad.set_xx(5, 0, w + dx / 2.0 - imprecision);
        pad.set_xx(5, 1, h / 2.0);
        pad.set_xx(6, 0, w / 2.0 + dx - imprecision);
        pad.set_xx(6, 1, h - imprecision);
        pad.set_xx(7, 0, 0.0 + dx / 2.0 + imprecision);
        pad.set_xx(7, 1, h / 2.0);
        assert!(pad.ok_xxt);
        if save_fig {
            pad.draw_shape_simple("/tmp/gemlab/gen_qua8.svg").unwrap();
        }
        pad
    }

    fn do_interpolate(pad: &mut Scratchpad, u_nodal: &Vector, gauss: &Gauss) -> Vector {
        let nnode = u_nodal.dim();
        let ngauss = gauss.npoint();
        let mut u_point = Vector::new(ngauss);
        for p in 0..ngauss {
            pad.calc_interp(gauss.coords(p));
            for m in 0..nnode {
                u_point[p] += pad.interp[m] * u_nodal[m];
            }
        }
        u_point
    }

    fn do_extrapolate(ee: &Matrix, u_point: &Vector) -> Vector {
        let nnode = ee.dims().0;
        let mut u_nodal = Vector::new(nnode);
        mat_vec_mul(&mut u_nodal, 1.0, &ee, &u_point).unwrap();
        u_nodal
    }

    #[test]
    fn get_extrap_matrix_works_qua4_ip4() {
        let mut pad = gen_qua4(20.0, 10.0, PI / 6.0, false);

        // extrapolation matrix
        let gauss = Gauss::new_sized(pad.kind.class(), 4).unwrap();
        let ee = get_extrap_matrix(&mut pad, &gauss).unwrap();

        // check
        let u_nodal_original = Vector::from(&[1.0, 2.0, 3.0, 4.0]); // original U values at nodes
        let u_point = do_interpolate(&mut pad, &u_nodal_original, &gauss); // interpolated U values @ integration points
        let u_nodal = do_extrapolate(&ee, &u_point);
        vec_approx_eq(&u_nodal, &u_nodal_original, 1e-15);
    }

    #[test]
    fn get_extrap_matrix_works_qua4_ip1() {
        let mut pad = gen_qua4(20.0, 10.0, PI / 6.0, false);

        // extrapolation matrix
        let gauss = Gauss::new_sized(pad.kind.class(), 1).unwrap();
        let ee = get_extrap_matrix(&mut pad, &gauss).unwrap();
        let ee_correct = Matrix::from(&[[1.0], [1.0], [1.0], [1.0]]); // the pseudo-inverse of P = [0.25, 0.25, 0.25, 0.25]
        mat_approx_eq(&ee, &ee_correct, 1e-15);

        // check
        let u_nodal_original = Vector::from(&[1.0, 2.0, 3.0, 4.0]); // original U values at nodes
        let u_point = do_interpolate(&mut pad, &u_nodal_original, &gauss); // interpolated U values @ integration points
        let u_nodal = do_extrapolate(&ee, &u_point);
        vec_approx_eq(&u_nodal, &u_nodal_original, 1.5); // we cannot get a better result with such low ngauss
    }

    #[test]
    fn get_extrap_matrix_works_qua8_ip9() {
        let mut pad = gen_qua8(20.0, 10.0, PI / 6.0, 1.0, false);

        // extrapolation matrix
        let gauss = Gauss::new_sized(pad.kind.class(), 9).unwrap();
        let ee = get_extrap_matrix(&mut pad, &gauss).unwrap();

        // check
        let u_nodal_original = Vector::from(&[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]); // original U values at nodes
        let u_point = do_interpolate(&mut pad, &u_nodal_original, &gauss); // interpolated U values @ integration points
        let u_nodal = do_extrapolate(&ee, &u_point);
        vec_approx_eq(&u_nodal, &u_nodal_original, 1e-13); // pretty good, with overdetermined system
    }

    #[test]
    fn get_extrap_matrix_works_qua8_ip4() {
        let mut pad = gen_qua8(20.0, 10.0, 0.0, 0.0, false);

        // extrapolation matrix
        let gauss = Gauss::new_sized(pad.kind.class(), 4).unwrap();
        let ee = get_extrap_matrix(&mut pad, &gauss).unwrap();

        // check
        let u_nodal_original = Vector::from(&[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]); // original U values at nodes
        let u_point = do_interpolate(&mut pad, &u_nodal_original, &gauss); // interpolated U values @ integration points
        let u_nodal = do_extrapolate(&ee, &u_point);
        vec_approx_eq(&u_nodal, &u_nodal_original, 12.3); // we cannot get better results with only 4 integ points
    }

    // Calculates E = P⁺
    fn calc_uncorrected_extrap_matrix(pad: &mut Scratchpad, gauss: &Gauss) -> Matrix {
        let nv = pad.deriv.dims().0;
        let np = gauss.npoint();
        let mut pp = get_interp_matrix(pad, gauss); // (np,nv)
        let mut pp_inv = Matrix::new(nv, np);
        mat_pseudo_inverse(&mut pp_inv, &mut pp).unwrap();
        return pp_inv;
    }

    #[test]
    fn get_extrap_matrix_works_durand_farias_example1() {
        let mut pad = gen_qua8(1.0, 1.0, 0.0, 0.0, false);

        // extrapolation matrix
        let gauss = Gauss::new_sized(pad.kind.class(), 4).unwrap();
        let ee = get_extrap_matrix(&mut pad, &gauss).unwrap();

        // extrapolate to nodes
        let u_point = Vector::from(&[0.5, 0.5, 0.5, 0.5]);
        let u_nodal = do_extrapolate(&ee, &u_point);

        // check
        let u_nodal_expected = Vector::from(&[0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]);
        vec_approx_eq(&u_nodal, &u_nodal_expected, 1e-14);

        if SAVE_FIGURE {
            let ee_uncorrected = calc_uncorrected_extrap_matrix(&mut pad, &gauss);
            let u_nodal_unc = do_extrapolate(&ee_uncorrected, &u_point);

            let n = 21;
            let (rr, ss, uu) = generate3d(-1.0, 1.0, -1.0, 1.0, n, n, |r, s| {
                pad.calc_interp(&[r, s]);
                vec_inner(&pad.interp, &u_nodal)
            });
            let (rr_unc, ss_unc, uu_unc) = generate3d(-1.0, 1.0, -1.0, 1.0, n, n, |r, s| {
                pad.calc_interp(&[r, s]);
                vec_inner(&pad.interp, &u_nodal_unc)
            });

            let mut surface = Surface::new();
            surface
                .set_with_surface(false)
                .set_with_wireframe(true)
                .set_wire_line_width(0.3)
                .draw(&rr, &ss, &uu);
            let mut surface_unc = Surface::new();

            surface_unc
                .set_with_surface(false)
                .set_with_wireframe(true)
                .set_wire_line_width(0.3)
                .draw(&rr_unc, &ss_unc, &uu_unc);

            let mut plot = Plot::new();
            plot.set_subplot_3d(1, 2, 1)
                .set_title("Corrected E matrix")
                .set_labels_3d("r", "s", "U")
                .set_zrange(-0.09, 0.7)
                .add(&surface)
                .set_subplot_3d(1, 2, 2)
                .set_title("Uncorrected E matrix")
                .set_labels_3d("r", "s", "U")
                .set_zrange(-0.09, 0.7)
                .add(&surface_unc)
                .set_save_pad_inches(0.25)
                .set_figure_size_points(600.0, 400.0)
                .save("/tmp/gemlab/test_durand_farias_example1.svg")
                .unwrap();
        }
    }

    #[test]
    fn get_extrap_matrix_works_durand_farias_example2() {
        let mut pad = gen_qua8(1.0, 1.0, 0.0, 0.0, false);

        // extrapolation matrix
        let gauss = Gauss::new_sized(pad.kind.class(), 4).unwrap();
        let ee = get_extrap_matrix(&mut pad, &gauss).unwrap();

        // set plane function at integration points
        let np = gauss.npoint();
        let x_ips = get_points_coords(&mut pad, &gauss).unwrap();
        let mut u_point = Vector::new(np);
        for p in 0..np {
            u_point[p] = x_ips[p][0] + x_ips[p][1] - 1.0;
        }

        // extrapolate to nodes
        let u_nodal = do_extrapolate(&ee, &u_point);

        // check
        for m in 0..u_nodal.dim() {
            let x = pad.xxt.get(0, m);
            let y = pad.xxt.get(1, m);
            approx_eq(x + y - 1.0, u_nodal[m], 1e-14);
        }

        if SAVE_FIGURE {
            let ee_uncorrected = calc_uncorrected_extrap_matrix(&mut pad, &gauss);
            let u_nodal_unc = do_extrapolate(&ee_uncorrected, &u_point);

            let n = 21;
            let (rr, ss, uu) = generate3d(-1.0, 1.0, -1.0, 1.0, n, n, |r, s| {
                pad.calc_interp(&[r, s]);
                vec_inner(&pad.interp, &u_nodal)
            });
            let (rr_unc, ss_unc, uu_unc) = generate3d(-1.0, 1.0, -1.0, 1.0, n, n, |r, s| {
                pad.calc_interp(&[r, s]);
                vec_inner(&pad.interp, &u_nodal_unc)
            });

            let mut surface = Surface::new();
            surface
                .set_with_surface(false)
                .set_with_wireframe(true)
                .set_wire_line_width(0.3)
                .draw(&rr, &ss, &uu);
            let mut surface_unc = Surface::new();

            surface_unc
                .set_with_surface(false)
                .set_with_wireframe(true)
                .set_wire_line_width(0.3)
                .draw(&rr_unc, &ss_unc, &uu_unc);

            let mut plot = Plot::new();
            plot.set_subplot_3d(1, 2, 1)
                .set_title("Corrected E matrix")
                .set_labels_3d("r", "s", "U")
                .set_zrange(-1.0, 1.0)
                .add(&surface)
                .set_subplot_3d(1, 2, 2)
                .set_title("Uncorrected E matrix")
                .set_labels_3d("r", "s", "U")
                .set_zrange(-1.0, 1.0)
                .add(&surface_unc)
                .set_save_pad_inches(0.25)
                .set_figure_size_points(600.0, 400.0)
                .save("/tmp/gemlab/test_durand_farias_example2.svg")
                .unwrap();
        }
    }

    fn check_hyperplane(pad: &mut Scratchpad, gauss: &Gauss) {
        // println!("{:?}: {} points", pad.kind, gauss.npoint());

        // set (hyper) plane function at integration points
        let geo_ndim = pad.kind.ndim();
        let np = gauss.npoint();
        let x_ips = get_points_coords(pad, &gauss).unwrap();
        let mut u_point = Vector::new(np);
        for p in 0..np {
            u_point[p] = -1.0;
            for d in 0..geo_ndim {
                u_point[p] += x_ips[p][d]; // -1 + x + y + z
            }
        }

        // perform extrapolation to nodes
        let ee = get_extrap_matrix(pad, &gauss).unwrap();
        let u_nodal = do_extrapolate(&ee, &u_point);

        // check
        for m in 0..u_nodal.dim() {
            let mut u_expected = -1.0;
            for d in 0..geo_ndim {
                u_expected += pad.xxt.get(d, m); // -1 + x + y + z
            }
            approx_eq(u_expected, u_nodal[m], 1e-13);
        }

        // interpolate back to integration points
        let mut u_point_interp = Vector::new(np);
        let pp = get_interp_matrix(pad, &gauss);
        mat_vec_mul(&mut u_point_interp, 1.0, &pp, &u_nodal).unwrap();

        // check
        vec_approx_eq(&u_point_interp, &u_point, 1e-14);
    }

    #[test]
    fn get_extrap_matrix_works_many_combos() {
        // for kind in [GeoKind::Tri6] {
        for kind in GeoKind::VALUES {
            // pad
            let space_ndim = usize::max(2, kind.ndim());
            let mut pad = Scratchpad::new(space_ndim, kind).unwrap();

            // set the real coordinates of nodes equal to their natural coordinates
            let (nnode, geo_ndim) = pad.deriv.dims();
            for m in 0..nnode {
                let r = kind.reference_coords(m);
                for j in 0..geo_ndim {
                    pad.set_xx(m, j, r[j]);
                }
                if geo_ndim == 1 {
                    pad.set_xx(m, 1, 123.456);
                }
            }

            // check with default integration points
            let gauss = Gauss::new(kind);
            check_hyperplane(&mut pad, &gauss);

            // check with minimum number of integration points
            let ngauss = match kind.class() {
                GeoClass::Lin => 2,
                GeoClass::Tri => 3,
                GeoClass::Qua => 4,
                GeoClass::Tet => 4,
                GeoClass::Hex => 8,
            };
            if ngauss != gauss.npoint() {
                let gauss_min = Gauss::new_sized(kind.class(), ngauss).unwrap();
                check_hyperplane(&mut pad, &gauss_min);
            }

            // check with maximum number of integration points
            let ngauss = match kind.class() {
                GeoClass::Lin => 5,
                GeoClass::Tri => 16,
                GeoClass::Qua => 16,
                GeoClass::Tet => 24,
                GeoClass::Hex => 64,
            };
            if ngauss != gauss.npoint() {
                let gauss_min = Gauss::new_sized(kind.class(), ngauss).unwrap();
                check_hyperplane(&mut pad, &gauss_min);
            }
        }
    }
}