gemlab 2.0.0

use crate::shapes::{GeoClass, GeoKind};
use crate::StrError;
use russell_lab::Matrix;

/// Holds a set of integration points ("Gauss") associated with a GeoClass
///
/// # 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>
pub struct Gauss {
    /// Holds the geometry class
    pub class: GeoClass,

    /// Holds the reference (natural) coordinates and weights in the array
    /// [[f64; 4]; NIP], where `NIP` is the number of integration points and
    /// `4` corresponds to the `r`, `s`, `t` coordinates and the weight `w`.
    ///
    /// In 1D, `r` and `t` are equal to zero. In 2D, `t` is equal to zero.
    ///
    /// The 1D and 2D geometries have `s` and/or `t` equal to zero, as appropriate.
    data: &'static [[f64; 4]],

    /// Holds the augmented natural coordinates matrix (for the extrapolation algorithm)
    ///
    /// See paragraph after Eq (30) of Reference # 1.
    pub ksi_hat: Option<Matrix>,

    /// Holds the pseudo-inverse matrix of the ksi_hat matrix (for the extrapolation algorithm)
    ///
    /// See paragraph after Eq (30) of Reference # 1.
    pub ksi_hat_inv: Option<Matrix>,
}

impl Gauss {
    /// Allocates a new instance with default number of integration points
    ///
    /// **Note:** The default set is decided based on the number of nodes in GeoKind
    /// and the possible use case in finite element analyses.
    ///
    /// # Examples
    ///
    /// ```
    /// use gemlab::integ::Gauss;
    /// use gemlab::shapes::GeoKind;
    ///
    /// let gauss = Gauss::new(GeoKind::Tet4);
    /// assert_eq!(gauss.npoint(), 4);
    /// ```
    pub fn new(kind: GeoKind) -> Self {
        let data: &'static [[f64; 4]] = match kind {
            // Lin
            GeoKind::Lin2 => &IP_LIN_LEGENDRE_2,
            GeoKind::Lin3 => &IP_LIN_LEGENDRE_3,
            GeoKind::Lin4 => &IP_LIN_LEGENDRE_4,
            GeoKind::Lin5 => &IP_LIN_LEGENDRE_5,
            // Tri
            GeoKind::Tri3 => &IP_TRI_INTERNAL_3,
            GeoKind::Tri6 => &IP_TRI_FELIPPA_7,
            GeoKind::Tri10 => &IP_TRI_INTERNAL_12,
            GeoKind::Tri15 => &IP_TRI_INTERNAL_16,
            // Qua
            GeoKind::Qua4 => &IP_QUA_LEGENDRE_4,
            GeoKind::Qua8 => &IP_QUA_LEGENDRE_9,
            GeoKind::Qua9 => &IP_QUA_LEGENDRE_9,
            GeoKind::Qua12 => &IP_QUA_LEGENDRE_16,
            GeoKind::Qua16 => &IP_QUA_LEGENDRE_16,
            GeoKind::Qua17 => &IP_QUA_LEGENDRE_16,
            // Tet
            GeoKind::Tet4 => &IP_TET_INTERNAL_4,
            GeoKind::Tet10 => &IP_TET_FELIPPA_14,
            GeoKind::Tet20 => &IP_TET_FELIPPA_24,
            // Hex
            GeoKind::Hex8 => &IP_HEX_LEGENDRE_8,
            GeoKind::Hex20 => &IP_HEX_LEGENDRE_27,
            GeoKind::Hex32 => &IP_HEX_LEGENDRE_64,
        };
        Gauss {
            class: kind.class(),
            data,
            ksi_hat: None,
            ksi_hat_inv: None,
        }
    }

    /// Allocates a new instance given the GeoClass and number of integration points
    ///
    /// # Input
    ///
    /// * `class` -- The geometry class
    /// * `ngauss` -- Number of integration points (see Options below)
    ///
    /// # Options
    ///
    /// `ngauss` for **Lin** class:
    ///
    /// * `1` -- Conventional Legendre integration points and weights
    /// * `2` -- Conventional Legendre integration points and weights
    /// * `3` -- Conventional Legendre integration points and weights
    /// * `4` -- Conventional Legendre integration points and weights
    /// * `5` -- Conventional Legendre integration points and weights
    ///
    /// `ngauss` for **Tri** class:
    ///
    /// * `1` -- Internal integration points and weights
    /// * `3` -- Internal integration points and weights
    /// * `4` -- Internal integration points and weights
    /// * `6` -- Internal integration points and weights (based on Felippa's code)
    /// * `7` -- Internal integration points and weights (based on Felippa's code)
    /// * `12` -- Internal integration points and weights
    /// * `16` -- Internal integration points and weights
    ///
    /// `ngauss` for **Qua** class:
    ///
    /// * `1` -- Conventional Legendre integration points and weights
    /// * `4` -- Conventional Legendre integration points and weights
    /// * `9` -- Conventional Legendre integration points and weights
    /// * `16` -- Conventional Legendre integration points and weights
    ///
    /// `ngauss` for **Tet** class:
    ///
    /// * `1` -- Internal integration points and weights, degree 1
    /// * `4` -- Internal integration points and weights, degree 2 (based on Felippa's code)
    /// * `5` -- Internal integration points and weights, degree 3 (with negative weight)
    /// * `8` -- Felippa's ips and weights, degree 3 (based on Felippa's code)
    /// * `14` -- Felippa's ips and weights, degree 4 (based on Felippa's code)
    /// * `15` -- Felippa's ips and weights, degree 5 (based on Felippa's code)
    /// * `24` -- Felippa's ips and weights, degree 6 (based on Felippa's code)
    ///
    /// `ngauss` for **Hex** class:
    ///
    /// * `6` -- Iron's integration points and weights
    /// * `8` -- Conventional Legendre integration points and weights, degree 3
    /// * `14` -- Iron's integration points and weights
    /// * `27` -- Conventional Legendre integration points and weights, degree 5
    /// * `64` -- Conventional Legendre integration points and weights, degree 7
    ///
    /// # Examples
    ///
    /// ```
    /// use gemlab::integ::Gauss;
    /// use gemlab::shapes::GeoClass;
    /// use gemlab::StrError;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     let gauss = Gauss::new_sized(GeoClass::Tet, 1)?;
    ///     assert_eq!(gauss.coords(0), &[0.25, 0.25, 0.25, 1.0/6.0]);
    ///     assert_eq!(gauss.weight(0), 1.0/6.0);
    ///     Ok(())
    /// }
    /// ```
    pub fn new_sized(class: GeoClass, ngauss: usize) -> Result<Self, StrError> {
        let data: &'static [[f64; 4]] = match class {
            // Lin
            GeoClass::Lin => match ngauss {
                1 => &IP_LIN_LEGENDRE_1,
                2 => &IP_LIN_LEGENDRE_2,
                3 => &IP_LIN_LEGENDRE_3,
                4 => &IP_LIN_LEGENDRE_4,
                5 => &IP_LIN_LEGENDRE_5,
                _ => return Err("requested number of integration points is not available for Lin class"),
            },
            // Tri
            GeoClass::Tri => match ngauss {
                1 => &IP_TRI_INTERNAL_1,
                3 => &IP_TRI_INTERNAL_3,
                4 => &IP_TRI_INTERNAL_4,
                6 => &IP_TRI_FELIPPA_6,
                7 => &IP_TRI_FELIPPA_7,
                12 => &IP_TRI_INTERNAL_12,
                16 => &IP_TRI_INTERNAL_16,
                _ => return Err("requested number of integration points is not available for Tri class"),
            },
            // Qua
            GeoClass::Qua => match ngauss {
                1 => &IP_QUA_LEGENDRE_1,
                4 => &IP_QUA_LEGENDRE_4,
                9 => &IP_QUA_LEGENDRE_9,
                16 => &IP_QUA_LEGENDRE_16,
                _ => return Err("requested number of integration points is not available for Qua class"),
            },
            // Tet
            GeoClass::Tet => match ngauss {
                1 => &IP_TET_INTERNAL_1,
                4 => &IP_TET_INTERNAL_4,
                5 => &IP_TET_INTERNAL_5,
                8 => &IP_TET_FELIPPA_8,
                14 => &IP_TET_FELIPPA_14,
                15 => &IP_TET_FELIPPA_15,
                24 => &IP_TET_FELIPPA_24,
                _ => return Err("requested number of integration points is not available for Tet class"),
            },
            // Hex
            GeoClass::Hex => match ngauss {
                6 => &IP_HEX_IRONS_6,
                8 => &IP_HEX_LEGENDRE_8,
                14 => &IP_HEX_IRONS_14,
                27 => &IP_HEX_LEGENDRE_27,
                64 => &IP_HEX_LEGENDRE_64,
                _ => return Err("requested number of integration points is not available for Hex class"),
            },
        };
        Ok(Gauss {
            class,
            data,
            ksi_hat: None,
            ksi_hat_inv: None,
        })
    }

    /// Allocates a new instance using either default or specified number of points
    ///
    /// # Input
    ///
    /// * `ngauss` -- Number of integration points (see options in [Gauss::new_sized()]).
    ///   If None, the default number is selected as described in [Gauss::new()]
    pub fn new_or_sized(kind: GeoKind, ngauss: Option<usize>) -> Result<Self, StrError> {
        match ngauss {
            Some(n) => Gauss::new_sized(kind.class(), n),
            None => Ok(Gauss::new(kind)),
        }
    }

    /// Returns the number of integration points
    #[inline]
    pub fn npoint(&self) -> usize {
        self.data.len()
    }

    /// Returns an access to the natural (reference) coordinates
    ///
    /// # Input
    ///
    /// * `p` -- the index of the integration point; must be in `0 ≤ p < npoint`
    ///   Otherwise a panic will occur.
    ///
    /// # Output
    ///
    /// * `coords` -- the `r`, `s`, `t` coordinates (followed by the weight---that can be ignored)
    ///   Note that, in 1D, `r` and `t` are equal to zero, and in 2D, `t` is equal to zero.
    ///
    /// # Panics
    ///
    /// A panic will occur if `p` is not in `0 ≤ p < npoint`
    #[inline]
    pub fn coords(&self, p: usize) -> &[f64; 4] {
        &self.data[p]
    }

    /// Returns the weight
    ///
    /// # Input
    ///
    /// * `p` -- the index of the integration point; must be in `0 ≤ p < npoint`
    ///   Otherwise a panic will occur.
    ///
    /// # Output
    ///
    /// * `weight` -- the weight associated with the integration point
    ///
    /// # Panics
    ///
    /// A panic will occur if `p` is not in `0 ≤ p < npoint`
    #[inline]
    pub fn weight(&self, p: usize) -> f64 {
        self.data[p][3]
    }
}

//
// Next, we define several constants with integration points data.
//

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

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_LIN_LEGENDRE_1: [[f64; 4]; 1] = [
    [0.0, 0.0, 0.0, 2.0]
];

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_LIN_LEGENDRE_2: [[f64; 4]; 2] = [
    [-0.5773502691896257, 0.0, 0.0, 1.0],
    [ 0.5773502691896257, 0.0, 0.0, 1.0],
];

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_LIN_LEGENDRE_3: [[f64; 4]; 3] = [
    [-0.7745966692414834, 0.0, 0.0, 0.5555555555555556],
    [ 0.0000000000000000, 0.0, 0.0, 0.8888888888888888],
    [ 0.7745966692414834, 0.0, 0.0, 0.5555555555555556],
];

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_LIN_LEGENDRE_4: [[f64; 4]; 4] = [
    [-0.8611363115940526, 0.0, 0.0, 0.3478548451374538],
    [-0.3399810435848562, 0.0, 0.0, 0.6521451548625462],
    [ 0.3399810435848562, 0.0, 0.0, 0.6521451548625462],
    [ 0.8611363115940526, 0.0, 0.0, 0.3478548451374538],
];

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_LIN_LEGENDRE_5: [[f64; 4]; 5] = [
    [-0.9061798459386640, 0.0, 0.0, 0.2369268850561891],
    [-0.5384693101056831, 0.0, 0.0, 0.4786286704993665],
    [ 0.0000000000000000, 0.0, 0.0, 0.5688888888888889],
    [ 0.5384693101056831, 0.0, 0.0, 0.4786286704993665],
    [ 0.9061798459386640, 0.0, 0.0, 0.2369268850561891],
];

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

/// Internal integration points and weights
#[rustfmt::skip]
const IP_TRI_INTERNAL_1: [[f64; 4]; 1] = [
    [1.0 / 3.0, 1.0 / 3.0, 0.0, 1.0 / 2.0],
];

/// Internal integration points and weights
#[rustfmt::skip]
const IP_TRI_INTERNAL_3: [[f64; 4]; 3] = [
    [1.0 / 6.0, 1.0 / 6.0, 0.0, 1.0 / 6.0],
    [2.0 / 3.0, 1.0 / 6.0, 0.0, 1.0 / 6.0],
    [1.0 / 6.0, 2.0 / 3.0, 0.0, 1.0 / 6.0],
];

/// Internal integration points and weights
#[rustfmt::skip]
const IP_TRI_INTERNAL_4: [[f64; 4]; 4] = [
    [1.0 / 3.0, 1.0 / 3.0, 0.0, -27.0 / 96.0],
    [1.0 / 5.0, 1.0 / 5.0, 0.0,  25.0 / 96.0],
    [3.0 / 5.0, 1.0 / 5.0, 0.0,  25.0 / 96.0],
    [1.0 / 5.0, 3.0 / 5.0, 0.0,  25.0 / 96.0],
];

/// Internal integration points and weights, 6 points
#[rustfmt::skip]
const IP_TRI_FELIPPA_6: [[f64; 4]; 6] = [
    [0.44594849091596488632,  0.44594849091596488632,  0.0, 0.11169079483900573285 ],
    [0.10810301816807022736,  0.44594849091596488632,  0.0, 0.11169079483900573285 ],
    [0.44594849091596488632,  0.10810301816807022736,  0.0, 0.11169079483900573285 ],
    [0.091576213509770743460, 0.091576213509770743460, 0.0, 0.054975871827660933819],
    [0.81684757298045851308,  0.091576213509770743460, 0.0, 0.054975871827660933819],
    [0.091576213509770743460, 0.81684757298045851308,  0.0, 0.054975871827660933819],
];

/// Internal integration points and weights, 7 points
#[rustfmt::skip]
const IP_TRI_FELIPPA_7: [[f64; 4]; 7] = [
    [0.10128650732345633880,  0.10128650732345633880,  0.0, 0.062969590272413576298],
    [0.79742698535308732240,  0.10128650732345633880,  0.0, 0.062969590272413576298],
    [0.10128650732345633880,  0.79742698535308732240,  0.0, 0.062969590272413576298],
    [0.47014206410511508977,  0.47014206410511508977,  0.0, 0.066197076394253090369],
    [0.059715871789769820459, 0.47014206410511508977,  0.0, 0.066197076394253090369],
    [0.47014206410511508977,  0.059715871789769820459, 0.0, 0.066197076394253090369],
    [0.33333333333333333333,  0.33333333333333333333,  0.0, 0.11250000000000000000 ],
];

/// Internal integration points and weights
#[rustfmt::skip]
const IP_TRI_INTERNAL_12: [[f64; 4]; 12] = [
    [0.873821971016996, 0.063089014491502, 0.0, 0.0254224531851035],
    [0.063089014491502, 0.873821971016996, 0.0, 0.0254224531851035],
    [0.063089014491502, 0.063089014491502, 0.0, 0.0254224531851035],
    [0.501426509658179, 0.249286745170910, 0.0, 0.0583931378631895],
    [0.249286745170910, 0.501426509658179, 0.0, 0.0583931378631895],
    [0.249286745170910, 0.249286745170910, 0.0, 0.0583931378631895],
    [0.053145049844817, 0.310352451033784, 0.0, 0.041425537809187 ],
    [0.310352451033784, 0.053145049844817, 0.0, 0.041425537809187 ],
    [0.053145049844817, 0.636502499121398, 0.0, 0.041425537809187 ],
    [0.310352451033784, 0.636502499121398, 0.0, 0.041425537809187 ],
    [0.636502499121398, 0.053145049844817, 0.0, 0.041425537809187 ],
    [0.636502499121398, 0.310352451033784, 0.0, 0.041425537809187 ],
];

/// Internal integration points and weights
#[rustfmt::skip]
const IP_TRI_INTERNAL_16: [[f64; 4]; 16] = [
    [3.33333333333333e-01, 3.33333333333333e-01, 0.0, 7.21578038388935e-02],
    [8.14148234145540e-02, 4.59292588292723e-01, 0.0, 4.75458171336425e-02],
    [4.59292588292723e-01, 8.14148234145540e-02, 0.0, 4.75458171336425e-02],
    [4.59292588292723e-01, 4.59292588292723e-01, 0.0, 4.75458171336425e-02],
    [6.58861384496480e-01, 1.70569307751760e-01, 0.0, 5.16086852673590e-02],
    [1.70569307751760e-01, 6.58861384496480e-01, 0.0, 5.16086852673590e-02],
    [1.70569307751760e-01, 1.70569307751760e-01, 0.0, 5.16086852673590e-02],
    [8.98905543365938e-01, 5.05472283170310e-02, 0.0, 1.62292488115990e-02],
    [5.05472283170310e-02, 8.98905543365938e-01, 0.0, 1.62292488115990e-02],
    [5.05472283170310e-02, 5.05472283170310e-02, 0.0, 1.62292488115990e-02],
    [8.39477740995800e-03, 2.63112829634638e-01, 0.0, 1.36151570872175e-02],
    [7.28492392955404e-01, 8.39477740995800e-03, 0.0, 1.36151570872175e-02],
    [2.63112829634638e-01, 7.28492392955404e-01, 0.0, 1.36151570872175e-02],
    [8.39477740995800e-03, 7.28492392955404e-01, 0.0, 1.36151570872175e-02],
    [7.28492392955404e-01, 2.63112829634638e-01, 0.0, 1.36151570872175e-02],
    [2.63112829634638e-01, 8.39477740995800e-03, 0.0, 1.36151570872175e-02],
];

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

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_QUA_LEGENDRE_1: [[f64; 4]; 1] = [
    [0.0, 0.0, 0.0, 4.0],
];

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_QUA_LEGENDRE_4: [[f64; 4]; 4] = [
    [-0.5773502691896257, -0.5773502691896257, 0.0, 1.0],
    [ 0.5773502691896257, -0.5773502691896257, 0.0, 1.0],
    [-0.5773502691896257,  0.5773502691896257, 0.0, 1.0],
    [ 0.5773502691896257,  0.5773502691896257, 0.0, 1.0],
];

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_QUA_LEGENDRE_9: [[f64; 4]; 9] = [
    [-0.7745966692414834, -0.7745966692414834, 0.0, 25.0 / 81.0],
    [ 0.0000000000000000, -0.7745966692414834, 0.0, 40.0 / 81.0],
    [ 0.7745966692414834, -0.7745966692414834, 0.0, 25.0 / 81.0],
    [-0.7745966692414834,  0.0000000000000000, 0.0, 40.0 / 81.0],
    [ 0.0000000000000000,  0.0000000000000000, 0.0, 64.0 / 81.0],
    [ 0.7745966692414834,  0.0000000000000000, 0.0, 40.0 / 81.0],
    [-0.7745966692414834,  0.7745966692414834, 0.0, 25.0 / 81.0],
    [ 0.0000000000000000,  0.7745966692414834, 0.0, 40.0 / 81.0],
    [ 0.7745966692414834,  0.7745966692414834, 0.0, 25.0 / 81.0],
];

/// Conventional Legendre integration points and weights
#[rustfmt::skip]
const IP_QUA_LEGENDRE_16: [[f64; 4]; 16] = [
    [-0.8611363115940526, -0.8611363115940526, 0.0, 0.1210029932856019],
    [-0.3399810435848563, -0.8611363115940526, 0.0, 0.2268518518518519],
    [ 0.3399810435848563, -0.8611363115940526, 0.0, 0.2268518518518519],
    [ 0.8611363115940526, -0.8611363115940526, 0.0, 0.1210029932856019],
    [-0.8611363115940526, -0.3399810435848563, 0.0, 0.2268518518518519],
    [-0.3399810435848563, -0.3399810435848563, 0.0, 0.4252933030106947],
    [ 0.3399810435848563, -0.3399810435848563, 0.0, 0.4252933030106947],
    [ 0.8611363115940526, -0.3399810435848563, 0.0, 0.2268518518518519],
    [-0.8611363115940526,  0.3399810435848563, 0.0, 0.2268518518518519],
    [-0.3399810435848563,  0.3399810435848563, 0.0, 0.4252933030106947],
    [ 0.3399810435848563,  0.3399810435848563, 0.0, 0.4252933030106947],
    [ 0.8611363115940526,  0.3399810435848563, 0.0, 0.2268518518518519],
    [-0.8611363115940526,  0.8611363115940526, 0.0, 0.1210029932856019],
    [-0.3399810435848563,  0.8611363115940526, 0.0, 0.2268518518518519],
    [ 0.3399810435848563,  0.8611363115940526, 0.0, 0.2268518518518519],
    [ 0.8611363115940526,  0.8611363115940526, 0.0, 0.1210029932856019],
];

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

/// Internal integration points and weights, 1 point, degree 1
/// 
/// Reference: Wriggers (2008), page 12, table 4.5,
/// Nonlinear Finite Element Methods, Springer.
#[rustfmt::skip]
const IP_TET_INTERNAL_1: [[f64; 4]; 1] = [
    [1.0 / 4.0, 1.0 / 4.0, 1.0 / 4.0, 1.0 / 6.0],
];

/// Internal integration points and weights, 4 points, degree 2
/// 
/// References: Felippa Advanced FEM, page 17-21, Table 17.1. See also Felippa (2004) A compendium of
/// FEM integration formulas for symbolic work Engineering Computations, 21, 7/8, pg 867
#[rustfmt::skip]
const IP_TET_INTERNAL_4: [[f64; 4]; 4] = [
    [0.13819660112501051518, 0.13819660112501051518, 0.13819660112501051518, 0.041666666666666666667],
    [0.58541019662496845446, 0.13819660112501051518, 0.13819660112501051518, 0.041666666666666666667], 
    [0.13819660112501051518, 0.58541019662496845446, 0.13819660112501051518, 0.041666666666666666667],
    [0.13819660112501051518, 0.13819660112501051518, 0.58541019662496845446, 0.041666666666666666667],
];

/// Internal integration points and weights, 4 points, degree 3
/// 
/// Reference: Wriggers (2008), page 12, table 4.5,
/// Nonlinear Finite Element Methods, Springer.
#[rustfmt::skip]
const IP_TET_INTERNAL_5: [[f64; 4]; 5] = [
    [1.0 / 4.0, 1.0 / 4.0, 1.0 / 4.0, -2.0 / 15.0],
    [1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0,  3.0 / 40.0],
    [1.0 / 6.0, 1.0 / 6.0, 1.0 / 2.0,  3.0 / 40.0],
    [1.0 / 6.0, 1.0 / 2.0, 1.0 / 6.0,  3.0 / 40.0],
    [1.0 / 2.0, 1.0 / 6.0, 1.0 / 6.0,  3.0 / 40.0],
];

/// Internal integration points and weights, 8 points, degree 3
/// 
/// References: Felippa Advanced FEM, page 17-21, Table 17.1. See also Felippa (2004) A compendium of
/// FEM integration formulas for symbolic work Engineering Computations, 21, 7/8, pg 867
#[rustfmt::skip]
const IP_TET_FELIPPA_8: [[f64; 4]; 8] = [
    [0.32805469671142664734,  0.32805469671142664734,  0.32805469671142664734,  0.023087994418643690387], 
    [0.015835909865720057993, 0.32805469671142664734,  0.32805469671142664734,  0.023087994418643690387], 
    [0.32805469671142664734,  0.015835909865720057993, 0.32805469671142664734,  0.023087994418643690387], 
    [0.32805469671142664734,  0.32805469671142664734,  0.015835909865720057993, 0.023087994418643690387], 
    [0.10695227393293068277,  0.10695227393293068277,  0.10695227393293068277,  0.018578672248022976279], 
    [0.67914317820120795168,  0.10695227393293068277,  0.10695227393293068277,  0.018578672248022976279], 
    [0.10695227393293068277,  0.67914317820120795168,  0.10695227393293068277,  0.018578672248022976279], 
    [0.10695227393293068277,  0.10695227393293068277,  0.67914317820120795168,  0.018578672248022976279],
];

/// Internal integration points and weights, 14 points, degree 4
/// 
/// References: Felippa Advanced FEM, page 17-21, Table 17.1. See also Felippa (2004) A compendium of
/// FEM integration formulas for symbolic work Engineering Computations, 21, 7/8, pg 867
#[rustfmt::skip]
const IP_TET_FELIPPA_14: [[f64; 4]; 14] = [
    [0.092735250310891226402, 0.092735250310891226402, 0.092735250310891226402, 0.012248840519393658257 ], 
    [0.72179424906732632079,  0.092735250310891226402, 0.092735250310891226402, 0.012248840519393658257 ], 
    [0.092735250310891226402, 0.72179424906732632079,  0.092735250310891226402, 0.012248840519393658257 ], 
    [0.092735250310891226402, 0.092735250310891226402, 0.72179424906732632079,  0.012248840519393658257 ], 
    [0.31088591926330060980,  0.31088591926330060980,  0.31088591926330060980,  0.018781320953002641800 ], 
    [0.067342242210098170608, 0.31088591926330060980,  0.31088591926330060980,  0.018781320953002641800 ], 
    [0.31088591926330060980,  0.067342242210098170608, 0.31088591926330060980,  0.018781320953002641800 ], 
    [0.31088591926330060980,  0.31088591926330060980,  0.067342242210098170608, 0.018781320953002641800 ], 
    [0.045503704125649649492, 0.45449629587435035051,  0.45449629587435035051,  0.0070910034628469110730], 
    [0.45449629587435035051,  0.045503704125649649492, 0.45449629587435035051,  0.0070910034628469110730], 
    [0.45449629587435035051,  0.45449629587435035051,  0.045503704125649649492, 0.0070910034628469110730], 
    [0.045503704125649649492, 0.045503704125649649492, 0.45449629587435035051,  0.0070910034628469110730], 
    [0.045503704125649649492, 0.45449629587435035051,  0.045503704125649649492, 0.0070910034628469110730], 
    [0.45449629587435035051,  0.045503704125649649492, 0.045503704125649649492, 0.0070910034628469110730],
];

/// Internal integration points and weights, 15 points, degree 5
/// 
/// References: Felippa Advanced FEM, page 17-21, Table 17.1. See also Felippa (2004) A compendium of
/// FEM integration formulas for symbolic work Engineering Computations, 21, 7/8, pg 867
#[rustfmt::skip]
const IP_TET_FELIPPA_15: [[f64; 4]; 15] = [
    [0.091971078052723032789, 0.091971078052723032789, 0.091971078052723032789, 0.011989513963169770002 ], 
    [0.72408676584183090163,  0.091971078052723032789, 0.091971078052723032789, 0.011989513963169770002 ], 
    [0.091971078052723032789, 0.72408676584183090163,  0.091971078052723032789, 0.011989513963169770002 ], 
    [0.091971078052723032789, 0.091971078052723032789, 0.72408676584183090163,  0.011989513963169770002 ], 
    [0.31979362782962990839,  0.31979362782962990839,  0.31979362782962990839,  0.011511367871045397547 ], 
    [0.040619116511110274837, 0.31979362782962990839,  0.31979362782962990839,  0.011511367871045397547 ], 
    [0.31979362782962990839,  0.040619116511110274837, 0.31979362782962990839,  0.011511367871045397547 ], 
    [0.31979362782962990839,  0.31979362782962990839,  0.040619116511110274837, 0.011511367871045397547 ], 
    [0.44364916731037084426,  0.056350832689629155741, 0.056350832689629155741, 0.0088183421516754850088], 
    [0.056350832689629155741, 0.44364916731037084426,  0.056350832689629155741, 0.0088183421516754850088], 
    [0.056350832689629155741, 0.056350832689629155741, 0.44364916731037084426,  0.0088183421516754850088], 
    [0.44364916731037084426,  0.44364916731037084426,  0.056350832689629155741, 0.0088183421516754850088], 
    [0.44364916731037084426,  0.056350832689629155741, 0.44364916731037084426,  0.0088183421516754850088], 
    [0.056350832689629155741, 0.44364916731037084426,  0.44364916731037084426,  0.0088183421516754850088], 
    [0.25000000000000000000,  0.25000000000000000000,  0.25000000000000000000,  0.019753086419753086420 ],
];

/// Internal integration points and weights, 24 points, degree 6
/// 
/// References: Felippa Advanced FEM, page 17-21, Table 17.1. See also Felippa (2004) A compendium of
/// FEM integration formulas for symbolic work Engineering Computations, 21, 7/8, pg 867
#[rustfmt::skip]
const IP_TET_FELIPPA_24: [[f64; 4]; 24] = [
    [0.21460287125915202929,  0.21460287125915202929,  0.21460287125915202929,  0.0066537917096945820166], 
    [0.35619138622254391213,  0.21460287125915202929,  0.21460287125915202929,  0.0066537917096945820166], 
    [0.21460287125915202929,  0.35619138622254391213,  0.21460287125915202929,  0.0066537917096945820166], 
    [0.21460287125915202929,  0.21460287125915202929,  0.35619138622254391213,  0.0066537917096945820166], 
    [0.040673958534611353116, 0.040673958534611353116, 0.040673958534611353116, 0.0016795351758867738247], 
    [0.87797812439616594065,  0.040673958534611353116, 0.040673958534611353116, 0.0016795351758867738247], 
    [0.040673958534611353116, 0.87797812439616594065,  0.040673958534611353116, 0.0016795351758867738247], 
    [0.040673958534611353116, 0.040673958534611353116, 0.87797812439616594065,  0.0016795351758867738247], 
    [0.32233789014227551034,  0.32233789014227551034,  0.32233789014227551034,  0.0092261969239424536825], 
    [0.032986329573173468968, 0.32233789014227551034,  0.32233789014227551034,  0.0092261969239424536825], 
    [0.32233789014227551034,  0.032986329573173468968, 0.32233789014227551034,  0.0092261969239424536825], 
    [0.32233789014227551034,  0.32233789014227551034,  0.032986329573173468968, 0.0092261969239424536825], 
    [0.26967233145831580803,  0.063661001875017525299, 0.063661001875017525299, 0.0080357142857142857143], 
    [0.063661001875017525299, 0.26967233145831580803,  0.063661001875017525299, 0.0080357142857142857143], 
    [0.063661001875017525299, 0.063661001875017525299, 0.26967233145831580803,  0.0080357142857142857143], 
    [0.60300566479164914137,  0.26967233145831580803,  0.063661001875017525299, 0.0080357142857142857143], 
    [0.60300566479164914137,  0.063661001875017525299, 0.26967233145831580803,  0.0080357142857142857143], 
    [0.063661001875017525299, 0.60300566479164914137,  0.26967233145831580803,  0.0080357142857142857143], 
    [0.60300566479164914137,  0.063661001875017525299, 0.063661001875017525299, 0.0080357142857142857143], 
    [0.063661001875017525299, 0.60300566479164914137,  0.063661001875017525299, 0.0080357142857142857143], 
    [0.063661001875017525299, 0.063661001875017525299, 0.60300566479164914137,  0.0080357142857142857143], 
    [0.26967233145831580803,  0.60300566479164914137,  0.063661001875017525299, 0.0080357142857142857143], 
    [0.26967233145831580803,  0.063661001875017525299, 0.60300566479164914137,  0.0080357142857142857143], 
    [0.063661001875017525299, 0.26967233145831580803,  0.60300566479164914137,  0.0080357142857142857143],
];

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

/// Iron's integration points and weights
#[rustfmt::skip]
const IP_HEX_IRONS_6: [[f64; 4]; 6] = [
    [-1.0,  0.0,  0.0, 4.0 / 3.0],
    [ 1.0,  0.0,  0.0, 4.0 / 3.0],
    [ 0.0, -1.0,  0.0, 4.0 / 3.0],
    [ 0.0,  1.0,  0.0, 4.0 / 3.0],
    [ 0.0,  0.0, -1.0, 4.0 / 3.0],
    [ 0.0,  0.0,  1.0, 4.0 / 3.0],
];

/// Conventional Legendre integration points and weights, 8 points, degree 3
#[rustfmt::skip]
const IP_HEX_LEGENDRE_8: [[f64; 4]; 8] = [
    [-0.5773502691896257, -0.5773502691896257, -0.5773502691896257, 1.0],
    [ 0.5773502691896257, -0.5773502691896257, -0.5773502691896257, 1.0],
    [-0.5773502691896257,  0.5773502691896257, -0.5773502691896257, 1.0],
    [ 0.5773502691896257,  0.5773502691896257, -0.5773502691896257, 1.0],
    [-0.5773502691896257, -0.5773502691896257,  0.5773502691896257, 1.0],
    [ 0.5773502691896257, -0.5773502691896257,  0.5773502691896257, 1.0],
    [-0.5773502691896257,  0.5773502691896257,  0.5773502691896257, 1.0],
    [ 0.5773502691896257,  0.5773502691896257,  0.5773502691896257, 1.0],
];

/// Iron's integration points and weights, 14 points
#[rustfmt::skip]
const IP_HEX_IRONS_14: [[f64; 4]; 14] = [
    [ 0.7958224257542215,  0.0000000000000000,  0.0000000000000000,  0.8864265927977839],
    [-0.7958224257542215,  0.0000000000000000,  0.0000000000000000,  0.8864265927977839],
    [ 0.0000000000000000,  0.7958224257542215,  0.0000000000000000,  0.8864265927977839],
    [ 0.0000000000000000, -0.7958224257542215,  0.0000000000000000,  0.8864265927977839],
    [ 0.0000000000000000,  0.0000000000000000,  0.7958224257542215,  0.8864265927977839],
    [ 0.0000000000000000,  0.0000000000000000, -0.7958224257542215,  0.8864265927977839],
    [ 0.7587869106393281,  0.7587869106393281,  0.7587869106393281,  0.3351800554016621],
    [-0.7587869106393281,  0.7587869106393281,  0.7587869106393281,  0.3351800554016621],
    [ 0.7587869106393281, -0.7587869106393281,  0.7587869106393281,  0.3351800554016621],
    [-0.7587869106393281, -0.7587869106393281,  0.7587869106393281,  0.3351800554016621],
    [ 0.7587869106393281,  0.7587869106393281, -0.7587869106393281,  0.3351800554016621],
    [-0.7587869106393281,  0.7587869106393281, -0.7587869106393281,  0.3351800554016621],
    [ 0.7587869106393281, -0.7587869106393281, -0.7587869106393281,  0.3351800554016621],
    [-0.7587869106393281, -0.7587869106393281, -0.7587869106393281,  0.3351800554016621],
];

/// Conventional Legendre integration points and weights, 27 points, degree 5
#[rustfmt::skip]
const IP_HEX_LEGENDRE_27: [[f64; 4]; 27] = [
    [-0.774596669241483, -0.774596669241483, -0.774596669241483, 0.171467764060357],
    [ 0.000000000000000, -0.774596669241483, -0.774596669241483, 0.274348422496571],
    [ 0.774596669241483, -0.774596669241483, -0.774596669241483, 0.171467764060357],
    [-0.774596669241483,  0.000000000000000, -0.774596669241483, 0.274348422496571],
    [ 0.000000000000000,  0.000000000000000, -0.774596669241483, 0.438957475994513],
    [ 0.774596669241483,  0.000000000000000, -0.774596669241483, 0.274348422496571],
    [-0.774596669241483,  0.774596669241483, -0.774596669241483, 0.171467764060357],
    [ 0.000000000000000,  0.774596669241483, -0.774596669241483, 0.274348422496571],
    [ 0.774596669241483,  0.774596669241483, -0.774596669241483, 0.171467764060357],
    [-0.774596669241483, -0.774596669241483,  0.000000000000000, 0.274348422496571],
    [ 0.000000000000000, -0.774596669241483,  0.000000000000000, 0.438957475994513],
    [ 0.774596669241483, -0.774596669241483,  0.000000000000000, 0.274348422496571],
    [-0.774596669241483,  0.000000000000000,  0.000000000000000, 0.438957475994513],
    [ 0.000000000000000,  0.000000000000000,  0.000000000000000, 0.702331961591221],
    [ 0.774596669241483,  0.000000000000000,  0.000000000000000, 0.438957475994513],
    [-0.774596669241483,  0.774596669241483,  0.000000000000000, 0.274348422496571],
    [ 0.000000000000000,  0.774596669241483,  0.000000000000000, 0.438957475994513],
    [ 0.774596669241483,  0.774596669241483,  0.000000000000000, 0.274348422496571],
    [-0.774596669241483, -0.774596669241483,  0.774596669241483, 0.171467764060357],
    [ 0.000000000000000, -0.774596669241483,  0.774596669241483, 0.274348422496571],
    [ 0.774596669241483, -0.774596669241483,  0.774596669241483, 0.171467764060357],
    [-0.774596669241483,  0.000000000000000,  0.774596669241483, 0.274348422496571],
    [ 0.000000000000000,  0.000000000000000,  0.774596669241483, 0.438957475994513],
    [ 0.774596669241483,  0.000000000000000,  0.774596669241483, 0.274348422496571],
    [-0.774596669241483,  0.774596669241483,  0.774596669241483, 0.171467764060357],
    [ 0.000000000000000,  0.774596669241483,  0.774596669241483, 0.274348422496571],
    [ 0.774596669241483,  0.774596669241483,  0.774596669241483, 0.171467764060357],
];

/// Conventional Legendre integration points and weights, 64 points, degree 7
#[rustfmt::skip]
const IP_HEX_LEGENDRE_64: [[f64; 4]; 64] = [
    [-0.86113631159405257522, -0.86113631159405257522, -0.86113631159405257522, 0.04209147749053145454],
    [-0.33998104358485626480, -0.86113631159405257522, -0.86113631159405257522, 0.07891151579507055098],
    [ 0.33998104358485626480, -0.86113631159405257522, -0.86113631159405257522, 0.07891151579507055098],
    [ 0.86113631159405257522, -0.86113631159405257522, -0.86113631159405257522, 0.04209147749053145454],
    [-0.86113631159405257522, -0.33998104358485626480, -0.86113631159405257522, 0.07891151579507055098],
    [-0.33998104358485626480, -0.33998104358485626480, -0.86113631159405257522, 0.14794033605678130087],
    [ 0.33998104358485626480, -0.33998104358485626480, -0.86113631159405257522, 0.14794033605678130087],
    [ 0.86113631159405257522, -0.33998104358485626480, -0.86113631159405257522, 0.07891151579507055098],
    [-0.86113631159405257522,  0.33998104358485626480, -0.86113631159405257522, 0.07891151579507055098],
    [-0.33998104358485626480,  0.33998104358485626480, -0.86113631159405257522, 0.14794033605678130087],
    [ 0.33998104358485626480,  0.33998104358485626480, -0.86113631159405257522, 0.14794033605678130087],
    [ 0.86113631159405257522,  0.33998104358485626480, -0.86113631159405257522, 0.07891151579507055098],
    [-0.86113631159405257522,  0.86113631159405257522, -0.86113631159405257522, 0.04209147749053145454],
    [-0.33998104358485626480,  0.86113631159405257522, -0.86113631159405257522, 0.07891151579507055098],
    [ 0.33998104358485626480,  0.86113631159405257522, -0.86113631159405257522, 0.07891151579507055098],
    [ 0.86113631159405257522,  0.86113631159405257522, -0.86113631159405257522, 0.04209147749053145454],
    [-0.86113631159405257522, -0.86113631159405257522, -0.33998104358485626480, 0.07891151579507055098],
    [-0.33998104358485626480, -0.86113631159405257522, -0.33998104358485626480, 0.14794033605678130087],
    [ 0.33998104358485626480, -0.86113631159405257522, -0.33998104358485626480, 0.14794033605678130087],
    [ 0.86113631159405257522, -0.86113631159405257522, -0.33998104358485626480, 0.07891151579507055098],
    [-0.86113631159405257522, -0.33998104358485626480, -0.33998104358485626480, 0.14794033605678130087],
    [-0.33998104358485626480, -0.33998104358485626480, -0.33998104358485626480, 0.27735296695391298990],
    [ 0.33998104358485626480, -0.33998104358485626480, -0.33998104358485626480, 0.27735296695391298990],
    [ 0.86113631159405257522, -0.33998104358485626480, -0.33998104358485626480, 0.14794033605678130087],
    [-0.86113631159405257522,  0.33998104358485626480, -0.33998104358485626480, 0.14794033605678130087],
    [-0.33998104358485626480,  0.33998104358485626480, -0.33998104358485626480, 0.27735296695391298990],
    [ 0.33998104358485626480,  0.33998104358485626480, -0.33998104358485626480, 0.27735296695391298990],
    [ 0.86113631159405257522,  0.33998104358485626480, -0.33998104358485626480, 0.14794033605678130087],
    [-0.86113631159405257522,  0.86113631159405257522, -0.33998104358485626480, 0.07891151579507055098],
    [-0.33998104358485626480,  0.86113631159405257522, -0.33998104358485626480, 0.14794033605678130087],
    [ 0.33998104358485626480,  0.86113631159405257522, -0.33998104358485626480, 0.14794033605678130087],
    [ 0.86113631159405257522,  0.86113631159405257522, -0.33998104358485626480, 0.07891151579507055098],
    [-0.86113631159405257522, -0.86113631159405257522,  0.33998104358485626480, 0.07891151579507055098],
    [-0.33998104358485626480, -0.86113631159405257522,  0.33998104358485626480, 0.14794033605678130087],
    [ 0.33998104358485626480, -0.86113631159405257522,  0.33998104358485626480, 0.14794033605678130087],
    [ 0.86113631159405257522, -0.86113631159405257522,  0.33998104358485626480, 0.07891151579507055098],
    [-0.86113631159405257522, -0.33998104358485626480,  0.33998104358485626480, 0.14794033605678130087],
    [-0.33998104358485626480, -0.33998104358485626480,  0.33998104358485626480, 0.27735296695391298990],
    [ 0.33998104358485626480, -0.33998104358485626480,  0.33998104358485626480, 0.27735296695391298990],
    [ 0.86113631159405257522, -0.33998104358485626480,  0.33998104358485626480, 0.14794033605678130087],
    [-0.86113631159405257522,  0.33998104358485626480,  0.33998104358485626480, 0.14794033605678130087],
    [-0.33998104358485626480,  0.33998104358485626480,  0.33998104358485626480, 0.27735296695391298990],
    [ 0.33998104358485626480,  0.33998104358485626480,  0.33998104358485626480, 0.27735296695391298990],
    [ 0.86113631159405257522,  0.33998104358485626480,  0.33998104358485626480, 0.14794033605678130087],
    [-0.86113631159405257522,  0.86113631159405257522,  0.33998104358485626480, 0.07891151579507055098],
    [-0.33998104358485626480,  0.86113631159405257522,  0.33998104358485626480, 0.14794033605678130087],
    [ 0.33998104358485626480,  0.86113631159405257522,  0.33998104358485626480, 0.14794033605678130087],
    [ 0.86113631159405257522,  0.86113631159405257522,  0.33998104358485626480, 0.07891151579507055098],
    [-0.86113631159405257522, -0.86113631159405257522,  0.86113631159405257522, 0.04209147749053145454],
    [-0.33998104358485626480, -0.86113631159405257522,  0.86113631159405257522, 0.07891151579507055098],
    [ 0.33998104358485626480, -0.86113631159405257522,  0.86113631159405257522, 0.07891151579507055098],
    [ 0.86113631159405257522, -0.86113631159405257522,  0.86113631159405257522, 0.04209147749053145454],
    [-0.86113631159405257522, -0.33998104358485626480,  0.86113631159405257522, 0.07891151579507055098],
    [-0.33998104358485626480, -0.33998104358485626480,  0.86113631159405257522, 0.14794033605678130087],
    [ 0.33998104358485626480, -0.33998104358485626480,  0.86113631159405257522, 0.14794033605678130087],
    [ 0.86113631159405257522, -0.33998104358485626480,  0.86113631159405257522, 0.07891151579507055098],
    [-0.86113631159405257522,  0.33998104358485626480,  0.86113631159405257522, 0.07891151579507055098],
    [-0.33998104358485626480,  0.33998104358485626480,  0.86113631159405257522, 0.14794033605678130087],
    [ 0.33998104358485626480,  0.33998104358485626480,  0.86113631159405257522, 0.14794033605678130087],
    [ 0.86113631159405257522,  0.33998104358485626480,  0.86113631159405257522, 0.07891151579507055098],
    [-0.86113631159405257522,  0.86113631159405257522,  0.86113631159405257522, 0.04209147749053145454],
    [-0.33998104358485626480,  0.86113631159405257522,  0.86113631159405257522, 0.07891151579507055098],
    [ 0.33998104358485626480,  0.86113631159405257522,  0.86113631159405257522, 0.07891151579507055098],
    [ 0.86113631159405257522,  0.86113631159405257522,  0.86113631159405257522, 0.04209147749053145454],
];

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

#[cfg(test)]
mod tests {
    use super::Gauss;
    use crate::shapes::{GeoClass, GeoKind};

    #[test]
    fn new_works() {
        // Lin
        assert_eq!(Gauss::new(GeoKind::Lin2).data.len(), 2);
        assert_eq!(Gauss::new(GeoKind::Lin3).data.len(), 3);
        assert_eq!(Gauss::new(GeoKind::Lin4).data.len(), 4);
        assert_eq!(Gauss::new(GeoKind::Lin5).data.len(), 5);
        // Tri
        assert_eq!(Gauss::new(GeoKind::Tri3).data.len(), 3);
        assert_eq!(Gauss::new(GeoKind::Tri6).data.len(), 7);
        assert_eq!(Gauss::new(GeoKind::Tri10).data.len(), 12);
        assert_eq!(Gauss::new(GeoKind::Tri15).data.len(), 16);
        // Qua
        assert_eq!(Gauss::new(GeoKind::Qua4).data.len(), 4);
        assert_eq!(Gauss::new(GeoKind::Qua8).data.len(), 9);
        assert_eq!(Gauss::new(GeoKind::Qua9).data.len(), 9);
        assert_eq!(Gauss::new(GeoKind::Qua12).data.len(), 16);
        assert_eq!(Gauss::new(GeoKind::Qua16).data.len(), 16);
        assert_eq!(Gauss::new(GeoKind::Qua17).data.len(), 16);
        // Tet
        assert_eq!(Gauss::new(GeoKind::Tet4).data.len(), 4);
        assert_eq!(Gauss::new(GeoKind::Tet10).data.len(), 14);
        assert_eq!(Gauss::new(GeoKind::Tet20).data.len(), 24);
        // Hex
        assert_eq!(Gauss::new(GeoKind::Hex8).data.len(), 8);
        assert_eq!(Gauss::new(GeoKind::Hex20).data.len(), 27);
        assert_eq!(Gauss::new(GeoKind::Hex32).data.len(), 64);
    }

    #[test]
    fn new_sized_works() {
        // Lin
        for ngauss in [1, 2, 3, 4, 5] {
            let gauss = Gauss::new_sized(GeoClass::Lin, ngauss).unwrap();
            assert_eq!(gauss.class, GeoClass::Lin);
            assert_eq!(gauss.data.len(), ngauss);
        }
        assert_eq!(
            Gauss::new_sized(GeoClass::Lin, 100).err(),
            Some("requested number of integration points is not available for Lin class")
        );

        // Tri
        for ngauss in [1, 3, 4, 6, 7, 12, 16] {
            let gauss = Gauss::new_sized(GeoClass::Tri, ngauss).unwrap();
            assert_eq!(gauss.class, GeoClass::Tri);
            assert_eq!(gauss.data.len(), ngauss);
        }
        assert_eq!(
            Gauss::new_sized(GeoClass::Tri, 100).err(),
            Some("requested number of integration points is not available for Tri class")
        );

        // Qua
        for ngauss in [1, 4, 9, 16] {
            let gauss = Gauss::new_sized(GeoClass::Qua, ngauss).unwrap();
            assert_eq!(gauss.class, GeoClass::Qua);
            assert_eq!(gauss.data.len(), ngauss);
        }
        assert_eq!(
            Gauss::new_sized(GeoClass::Qua, 100).err(),
            Some("requested number of integration points is not available for Qua class")
        );

        // Tet
        for ngauss in [1, 4, 5, 8, 14, 15, 24] {
            let gauss = Gauss::new_sized(GeoClass::Tet, ngauss).unwrap();
            assert_eq!(gauss.class, GeoClass::Tet);
            assert_eq!(gauss.data.len(), ngauss);
        }
        assert_eq!(
            Gauss::new_sized(GeoClass::Tet, 100).err(),
            Some("requested number of integration points is not available for Tet class")
        );

        // Hex
        for ngauss in [6, 8, 14, 27, 64] {
            let gauss = Gauss::new_sized(GeoClass::Hex, ngauss).unwrap();
            assert_eq!(gauss.class, GeoClass::Hex);
            assert_eq!(gauss.data.len(), ngauss);
        }
        assert_eq!(
            Gauss::new_sized(GeoClass::Hex, 100).err(),
            Some("requested number of integration points is not available for Hex class")
        );
    }

    #[test]
    fn new_or_sized_works() {
        assert_eq!(Gauss::new_or_sized(GeoKind::Qua8, None).unwrap().npoint(), 9);
        assert_eq!(Gauss::new_or_sized(GeoKind::Qua8, Some(1)).unwrap().npoint(), 1);
    }
}