//! Finite Element Helpers for Structural Analysis
//!
//! This module provides element types, shape functions, and numerical
//! integration routines for finite element analysis.

use structural_shared::{Element2D, Element3D, Point2D, Point3D};

// ============================================================================
// Element Trait
// ============================================================================

/// Trait for finite element shape functions.
pub trait Element {
    /// Number of nodes in the element.
    fn num_nodes(&self) -> usize;

    /// Evaluate shape functions at natural coordinates.
    fn shape_functions(&self, xi: &[f64]) -> Vec<f64>;

    /// Evaluate shape function derivatives at natural coordinates.
    /// Returns (dN/d_xi, dN/d_eta, ...) for each node.
    fn shape_derivatives(&self, xi: &[f64]) -> Vec<Vec<f64>>;

    /// Get the number of spatial dimensions.
    fn dimensions(&self) -> usize;
}

// ============================================================================
// Triangle Element (3-node)
// ============================================================================

/// Linear triangle element (3 nodes).
#[derive(Debug, Clone, Copy)]
pub struct Triangle {
    /// Node coordinates.
    nodes: [Point2D; 3],
}

impl Triangle {
    /// Create a new triangle element.
    pub fn new(nodes: [Point2D; 3]) -> Self {
        Self { nodes }
    }

    /// Create from connectivity and vertex list.
    pub fn from_connectivity(connectivity: &[usize; 3], vertices: &[Point2D]) -> Self {
        Self {
            nodes: [
                vertices[connectivity[0]],
                vertices[connectivity[1]],
                vertices[connectivity[2]],
            ],
        }
    }

    /// Compute element area.
    pub fn area(&self) -> f64 {
        let [n0, n1, n2] = self.nodes;
        0.5 * ((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)).abs()
    }

    /// Compute element centroid.
    pub fn centroid(&self) -> Point2D {
        let [n0, n1, n2] = self.nodes;
        Point2D::new((n0.x + n1.x + n2.x) / 3.0, (n0.y + n1.y + n2.y) / 3.0)
    }

    /// Map natural coordinates to physical coordinates.
    pub fn map_to_physical(&self, xi: f64, eta: f64) -> Point2D {
        let n = self.shape_functions(&[xi, eta]);
        let [n0, n1, n2] = self.nodes;
        Point2D::new(
            n[0] * n0.x + n[1] * n1.x + n[2] * n2.x,
            n[0] * n0.y + n[1] * n1.y + n[2] * n2.y,
        )
    }

    /// Compute Jacobian matrix.
    pub fn jacobian(&self) -> [[f64; 2]; 2] {
        let [n0, n1, n2] = self.nodes;
        // For linear triangle, Jacobian is constant
        [[n1.x - n0.x, n2.x - n0.x], [n1.y - n0.y, n2.y - n0.y]]
    }

    /// Compute Jacobian determinant.
    pub fn jacobian_det(&self) -> f64 {
        let j = self.jacobian();
        j[0][0] * j[1][1] - j[0][1] * j[1][0]
    }

    /// Get nodes.
    pub fn nodes(&self) -> &[Point2D; 3] {
        &self.nodes
    }
}

impl Element for Triangle {
    fn num_nodes(&self) -> usize {
        3
    }

    fn shape_functions(&self, xi: &[f64]) -> Vec<f64> {
        // Natural coordinates for triangle: xi (L2), eta (L3)
        // L1 = 1 - xi - eta, L2 = xi, L3 = eta
        let (xi_val, eta_val) = (xi[0], xi[1]);
        vec![1.0 - xi_val - eta_val, xi_val, eta_val]
    }

    fn shape_derivatives(&self, _xi: &[f64]) -> Vec<Vec<f64>> {
        // dN/dxi, dN/deta for each node
        // N1 = 1 - xi - eta: dN1/dxi = -1, dN1/deta = -1
        // N2 = xi: dN2/dxi = 1, dN2/deta = 0
        // N3 = eta: dN3/dxi = 0, dN3/deta = 1
        vec![
            vec![-1.0, -1.0], // Node 1
            vec![1.0, 0.0],   // Node 2
            vec![0.0, 1.0],   // Node 3
        ]
    }

    fn dimensions(&self) -> usize {
        2
    }
}

// ============================================================================
// Quadrilateral Element (4-node)
// ============================================================================

/// Bilinear quadrilateral element (4 nodes).
#[derive(Debug, Clone, Copy)]
pub struct Quad {
    /// Node coordinates (counter-clockwise ordering).
    nodes: [Point2D; 4],
}

impl Quad {
    /// Create a new quad element.
    pub fn new(nodes: [Point2D; 4]) -> Self {
        Self { nodes }
    }

    /// Create from connectivity and vertex list.
    pub fn from_connectivity(connectivity: &[usize; 4], vertices: &[Point2D]) -> Self {
        Self {
            nodes: [
                vertices[connectivity[0]],
                vertices[connectivity[1]],
                vertices[connectivity[2]],
                vertices[connectivity[3]],
            ],
        }
    }

    /// Compute approximate element area.
    pub fn area(&self) -> f64 {
        let [n0, n1, n2, n3] = self.nodes;
        // Area = 0.5 * |diagonal1 x diagonal2|
        0.5 * ((n2.x - n0.x) * (n3.y - n1.y) - (n3.x - n1.x) * (n2.y - n0.y)).abs()
    }

    /// Compute element centroid.
    pub fn centroid(&self) -> Point2D {
        let [n0, n1, n2, n3] = self.nodes;
        Point2D::new(
            (n0.x + n1.x + n2.x + n3.x) / 4.0,
            (n0.y + n1.y + n2.y + n3.y) / 4.0,
        )
    }

    /// Map natural coordinates to physical coordinates.
    pub fn map_to_physical(&self, xi: f64, eta: f64) -> Point2D {
        let n = self.shape_functions(&[xi, eta]);
        let nodes = &self.nodes;
        Point2D::new(
            n[0] * nodes[0].x + n[1] * nodes[1].x + n[2] * nodes[2].x + n[3] * nodes[3].x,
            n[0] * nodes[0].y + n[1] * nodes[1].y + n[2] * nodes[2].y + n[3] * nodes[3].y,
        )
    }

    /// Compute Jacobian matrix at natural coordinates.
    pub fn jacobian(&self, xi: f64, eta: f64) -> [[f64; 2]; 2] {
        let dn = self.shape_derivatives(&[xi, eta]);
        let nodes = &self.nodes;

        let mut j = [[0.0; 2]; 2];
        for i in 0..4 {
            j[0][0] += dn[i][0] * nodes[i].x;
            j[0][1] += dn[i][0] * nodes[i].y;
            j[1][0] += dn[i][1] * nodes[i].x;
            j[1][1] += dn[i][1] * nodes[i].y;
        }
        j
    }

    /// Compute Jacobian determinant at natural coordinates.
    pub fn jacobian_det(&self, xi: f64, eta: f64) -> f64 {
        let j = self.jacobian(xi, eta);
        j[0][0] * j[1][1] - j[0][1] * j[1][0]
    }

    /// Get nodes.
    pub fn nodes(&self) -> &[Point2D; 4] {
        &self.nodes
    }
}

impl Element for Quad {
    fn num_nodes(&self) -> usize {
        4
    }

    fn shape_functions(&self, xi: &[f64]) -> Vec<f64> {
        // Bilinear shape functions on [-1, 1] x [-1, 1]
        let (xi_val, eta_val) = (xi[0], xi[1]);
        vec![
            0.25 * (1.0 - xi_val) * (1.0 - eta_val), // N1 at (-1, -1)
            0.25 * (1.0 + xi_val) * (1.0 - eta_val), // N2 at (1, -1)
            0.25 * (1.0 + xi_val) * (1.0 + eta_val), // N3 at (1, 1)
            0.25 * (1.0 - xi_val) * (1.0 + eta_val), // N4 at (-1, 1)
        ]
    }

    fn shape_derivatives(&self, xi: &[f64]) -> Vec<Vec<f64>> {
        let (xi_val, eta_val) = (xi[0], xi[1]);
        vec![
            vec![-0.25 * (1.0 - eta_val), -0.25 * (1.0 - xi_val)], // dN1
            vec![0.25 * (1.0 - eta_val), -0.25 * (1.0 + xi_val)],  // dN2
            vec![0.25 * (1.0 + eta_val), 0.25 * (1.0 + xi_val)],   // dN3
            vec![-0.25 * (1.0 + eta_val), 0.25 * (1.0 - xi_val)],  // dN4
        ]
    }

    fn dimensions(&self) -> usize {
        2
    }
}

// ============================================================================
// Tetrahedron Element (4-node)
// ============================================================================

/// Linear tetrahedron element (4 nodes).
#[derive(Debug, Clone, Copy)]
pub struct Tetrahedron {
    /// Node coordinates.
    nodes: [Point3D; 4],
}

impl Tetrahedron {
    /// Create a new tetrahedron element.
    pub fn new(nodes: [Point3D; 4]) -> Self {
        Self { nodes }
    }

    /// Create from connectivity and vertex list.
    pub fn from_connectivity(connectivity: &[usize; 4], vertices: &[Point3D]) -> Self {
        Self {
            nodes: [
                vertices[connectivity[0]],
                vertices[connectivity[1]],
                vertices[connectivity[2]],
                vertices[connectivity[3]],
            ],
        }
    }

    /// Compute element volume.
    pub fn volume(&self) -> f64 {
        let [n0, n1, n2, n3] = self.nodes;
        let v10 = (n1.x - n0.x, n1.y - n0.y, n1.z - n0.z);
        let v20 = (n2.x - n0.x, n2.y - n0.y, n2.z - n0.z);
        let v30 = (n3.x - n0.x, n3.y - n0.y, n3.z - n0.z);

        // Volume = |det([v10, v20, v30])| / 6
        let det = v10.0 * (v20.1 * v30.2 - v20.2 * v30.1) - v10.1 * (v20.0 * v30.2 - v20.2 * v30.0)
            + v10.2 * (v20.0 * v30.1 - v20.1 * v30.0);

        det.abs() / 6.0
    }

    /// Compute element centroid.
    pub fn centroid(&self) -> Point3D {
        let [n0, n1, n2, n3] = self.nodes;
        Point3D::new(
            (n0.x + n1.x + n2.x + n3.x) / 4.0,
            (n0.y + n1.y + n2.y + n3.y) / 4.0,
            (n0.z + n1.z + n2.z + n3.z) / 4.0,
        )
    }

    /// Get nodes.
    pub fn nodes(&self) -> &[Point3D; 4] {
        &self.nodes
    }
}

impl Element for Tetrahedron {
    fn num_nodes(&self) -> usize {
        4
    }

    fn shape_functions(&self, xi: &[f64]) -> Vec<f64> {
        // Natural coordinates: xi (L2), eta (L3), zeta (L4)
        // L1 = 1 - xi - eta - zeta
        let (xi_val, eta_val, zeta_val) = (xi[0], xi[1], xi[2]);
        vec![1.0 - xi_val - eta_val - zeta_val, xi_val, eta_val, zeta_val]
    }

    fn shape_derivatives(&self, _xi: &[f64]) -> Vec<Vec<f64>> {
        // dN/dxi, dN/deta, dN/dzeta for each node
        vec![
            vec![-1.0, -1.0, -1.0], // Node 1
            vec![1.0, 0.0, 0.0],    // Node 2
            vec![0.0, 1.0, 0.0],    // Node 3
            vec![0.0, 0.0, 1.0],    // Node 4
        ]
    }

    fn dimensions(&self) -> usize {
        3
    }
}

// ============================================================================
// Gauss Quadrature
// ============================================================================

/// Gauss quadrature points and weights.
#[derive(Debug, Clone)]
pub struct GaussQuadrature {
    /// Quadrature points (natural coordinates).
    pub points: Vec<Vec<f64>>,
    /// Quadrature weights.
    pub weights: Vec<f64>,
}

impl GaussQuadrature {
    /// Create 1-point quadrature for triangles.
    pub fn triangle_1() -> Self {
        Self {
            points: vec![vec![1.0 / 3.0, 1.0 / 3.0]],
            weights: vec![0.5], // Area of reference triangle
        }
    }

    /// Create 3-point quadrature for triangles.
    pub fn triangle_3() -> Self {
        Self {
            points: vec![
                vec![1.0 / 6.0, 1.0 / 6.0],
                vec![2.0 / 3.0, 1.0 / 6.0],
                vec![1.0 / 6.0, 2.0 / 3.0],
            ],
            weights: vec![1.0 / 6.0; 3],
        }
    }

    /// Create 4-point (2x2) quadrature for quads.
    pub fn quad_2x2() -> Self {
        let c = 1.0 / 3.0f64.sqrt();
        Self {
            points: vec![vec![-c, -c], vec![c, -c], vec![c, c], vec![-c, c]],
            weights: vec![1.0; 4],
        }
    }

    /// Create 9-point (3x3) quadrature for quads.
    pub fn quad_3x3() -> Self {
        let c = (3.0 / 5.0f64).sqrt();
        let w1 = 5.0 / 9.0;
        let w2 = 8.0 / 9.0;

        let mut points = Vec::with_capacity(9);
        let mut weights = Vec::with_capacity(9);

        let coords = [-c, 0.0, c];
        let wts = [w1, w2, w1];

        for (j, &eta) in coords.iter().enumerate() {
            for (i, &xi) in coords.iter().enumerate() {
                points.push(vec![xi, eta]);
                weights.push(wts[i] * wts[j]);
            }
        }

        Self { points, weights }
    }

    /// Create 1-point quadrature for tetrahedra.
    pub fn tetrahedron_1() -> Self {
        Self {
            points: vec![vec![0.25, 0.25, 0.25]],
            weights: vec![1.0 / 6.0], // Volume of reference tetrahedron
        }
    }

    /// Create 4-point quadrature for tetrahedra.
    pub fn tetrahedron_4() -> Self {
        let a = (5.0 - 5.0f64.sqrt()) / 20.0;
        let b = (5.0 + 3.0 * 5.0f64.sqrt()) / 20.0;

        Self {
            points: vec![vec![a, a, a], vec![b, a, a], vec![a, b, a], vec![a, a, b]],
            weights: vec![1.0 / 24.0; 4],
        }
    }

    /// Number of integration points.
    pub fn num_points(&self) -> usize {
        self.points.len()
    }

    /// Integrate a function over the element.
    pub fn integrate<F>(&self, func: F) -> f64
    where
        F: Fn(&[f64]) -> f64,
    {
        let mut result = 0.0;
        for (i, point) in self.points.iter().enumerate() {
            result += self.weights[i] * func(point);
        }
        result
    }
}

// ============================================================================
// B-Matrix (Strain-Displacement)
// ============================================================================

/// Compute B-matrix (strain-displacement) for 2D elements.
pub fn compute_b_matrix_2d(
    shape_derivatives: &[Vec<f64>],
    jacobian_inv: &[[f64; 2]; 2],
) -> Vec<Vec<f64>> {
    let num_nodes = shape_derivatives.len();
    // B-matrix: 3 rows (eps_xx, eps_yy, gamma_xy), 2*num_nodes columns

    let mut b = vec![vec![0.0; 2 * num_nodes]; 3];

    for i in 0..num_nodes {
        // Transform derivatives to physical coordinates
        let dn_dx = jacobian_inv[0][0] * shape_derivatives[i][0]
            + jacobian_inv[0][1] * shape_derivatives[i][1];
        let dn_dy = jacobian_inv[1][0] * shape_derivatives[i][0]
            + jacobian_inv[1][1] * shape_derivatives[i][1];

        // B-matrix entries
        // eps_xx = du/dx
        b[0][2 * i] = dn_dx;
        // eps_yy = dv/dy
        b[1][2 * i + 1] = dn_dy;
        // gamma_xy = du/dy + dv/dx
        b[2][2 * i] = dn_dy;
        b[2][2 * i + 1] = dn_dx;
    }

    b
}

/// Compute inverse of 2x2 matrix.
pub fn invert_2x2(m: &[[f64; 2]; 2]) -> Option<[[f64; 2]; 2]> {
    let det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
    if det.abs() < 1e-15 {
        return None;
    }
    let inv_det = 1.0 / det;
    Some([
        [m[1][1] * inv_det, -m[0][1] * inv_det],
        [-m[1][0] * inv_det, m[0][0] * inv_det],
    ])
}

// ============================================================================
// Element Factory
// ============================================================================

/// Create element from Element2D enum.
pub fn create_element_2d(elem: &Element2D, vertices: &[Point2D]) -> Box<dyn Element> {
    match elem {
        Element2D::Triangle(conn) => Box::new(Triangle::from_connectivity(conn, vertices)),
        Element2D::Quad(conn) => Box::new(Quad::from_connectivity(conn, vertices)),
        Element2D::Triangle6(_) => {
            // Simplified: use linear triangle
            Box::new(Triangle::new([vertices[0], vertices[1], vertices[2]]))
        }
        Element2D::Quad8(_) => {
            // Simplified: use bilinear quad
            Box::new(Quad::new([
                vertices[0],
                vertices[1],
                vertices[2],
                vertices[3],
            ]))
        }
    }
}

/// Create element from Element3D enum.
pub fn create_element_3d(elem: &Element3D, vertices: &[Point3D]) -> Box<dyn Element> {
    match elem {
        Element3D::Tetrahedron(conn) => Box::new(Tetrahedron::from_connectivity(conn, vertices)),
        _ => {
            // Simplified: use linear tetrahedron for other types
            Box::new(Tetrahedron::new([
                vertices[0],
                vertices[1],
                vertices[2],
                vertices[3],
            ]))
        }
    }
}

// ============================================================================
// Tests
// ============================================================================

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

    #[test]
    fn test_triangle_area() {
        let tri = Triangle::new([
            Point2D::new(0.0, 0.0),
            Point2D::new(1.0, 0.0),
            Point2D::new(0.0, 1.0),
        ]);
        assert!((tri.area() - 0.5).abs() < 1e-10);
    }

    #[test]
    fn test_triangle_centroid() {
        let tri = Triangle::new([
            Point2D::new(0.0, 0.0),
            Point2D::new(3.0, 0.0),
            Point2D::new(0.0, 3.0),
        ]);
        let c = tri.centroid();
        assert!((c.x - 1.0).abs() < 1e-10);
        assert!((c.y - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_triangle_shape_functions() {
        let tri = Triangle::new([
            Point2D::new(0.0, 0.0),
            Point2D::new(1.0, 0.0),
            Point2D::new(0.0, 1.0),
        ]);

        // At node 1 (0, 0): N1 = 1, N2 = 0, N3 = 0
        let n = tri.shape_functions(&[0.0, 0.0]);
        assert!((n[0] - 1.0).abs() < 1e-10);
        assert!(n[1].abs() < 1e-10);
        assert!(n[2].abs() < 1e-10);

        // At centroid: N1 = N2 = N3 = 1/3
        let n = tri.shape_functions(&[1.0 / 3.0, 1.0 / 3.0]);
        assert!((n[0] - 1.0 / 3.0).abs() < 1e-10);
        assert!((n[1] - 1.0 / 3.0).abs() < 1e-10);
        assert!((n[2] - 1.0 / 3.0).abs() < 1e-10);
    }

    #[test]
    fn test_quad_area() {
        let quad = Quad::new([
            Point2D::new(0.0, 0.0),
            Point2D::new(2.0, 0.0),
            Point2D::new(2.0, 1.0),
            Point2D::new(0.0, 1.0),
        ]);
        assert!((quad.area() - 2.0).abs() < 1e-10);
    }

    #[test]
    fn test_quad_shape_functions() {
        let quad = Quad::new([
            Point2D::new(0.0, 0.0),
            Point2D::new(1.0, 0.0),
            Point2D::new(1.0, 1.0),
            Point2D::new(0.0, 1.0),
        ]);

        // At corner (-1, -1): N1 = 1, others = 0
        let n = quad.shape_functions(&[-1.0, -1.0]);
        assert!((n[0] - 1.0).abs() < 1e-10);
        assert!(n[1].abs() < 1e-10);
        assert!(n[2].abs() < 1e-10);
        assert!(n[3].abs() < 1e-10);

        // At center (0, 0): all = 0.25
        let n = quad.shape_functions(&[0.0, 0.0]);
        for ni in &n {
            assert!((*ni - 0.25).abs() < 1e-10);
        }
    }

    #[test]
    fn test_tetrahedron_volume() {
        let tet = Tetrahedron::new([
            Point3D::new(0.0, 0.0, 0.0),
            Point3D::new(1.0, 0.0, 0.0),
            Point3D::new(0.0, 1.0, 0.0),
            Point3D::new(0.0, 0.0, 1.0),
        ]);
        assert!((tet.volume() - 1.0 / 6.0).abs() < 1e-10);
    }

    #[test]
    fn test_tetrahedron_shape_functions() {
        let tet = Tetrahedron::new([
            Point3D::new(0.0, 0.0, 0.0),
            Point3D::new(1.0, 0.0, 0.0),
            Point3D::new(0.0, 1.0, 0.0),
            Point3D::new(0.0, 0.0, 1.0),
        ]);

        // At node 1 (0, 0, 0): N1 = 1, others = 0
        let n = tet.shape_functions(&[0.0, 0.0, 0.0]);
        assert!((n[0] - 1.0).abs() < 1e-10);
        assert!(n[1].abs() < 1e-10);
        assert!(n[2].abs() < 1e-10);
        assert!(n[3].abs() < 1e-10);
    }

    #[test]
    fn test_gauss_triangle_1() {
        let quad = GaussQuadrature::triangle_1();
        assert_eq!(quad.num_points(), 1);
        assert!((quad.weights[0] - 0.5).abs() < 1e-10);
    }

    #[test]
    fn test_gauss_triangle_3() {
        let quad = GaussQuadrature::triangle_3();
        assert_eq!(quad.num_points(), 3);
        let sum: f64 = quad.weights.iter().sum();
        assert!((sum - 0.5).abs() < 1e-10); // Sum should equal reference area
    }

    #[test]
    fn test_gauss_quad_2x2() {
        let quad = GaussQuadrature::quad_2x2();
        assert_eq!(quad.num_points(), 4);
        let sum: f64 = quad.weights.iter().sum();
        assert!((sum - 4.0).abs() < 1e-10); // Sum should equal reference area (2x2)
    }

    #[test]
    fn test_gauss_integrate() {
        let quad = GaussQuadrature::quad_2x2();
        // Integrate f(x,y) = 1 over [-1,1]x[-1,1], should give 4
        let result = quad.integrate(|_| 1.0);
        assert!((result - 4.0).abs() < 1e-10);
    }

    #[test]
    fn test_gauss_integrate_polynomial() {
        let quad = GaussQuadrature::quad_2x2();
        // Integrate f(x,y) = x^2 over [-1,1]x[-1,1]
        // Exact: integral of x^2 from -1 to 1 = 2/3, times 2 (y range) = 4/3
        let result = quad.integrate(|xi| xi[0] * xi[0]);
        assert!((result - 4.0 / 3.0).abs() < 0.1); // 2x2 is exact for degree 3
    }

    #[test]
    fn test_invert_2x2() {
        let m = [[2.0, 1.0], [1.0, 2.0]];
        let inv = invert_2x2(&m).unwrap();
        // Check M * M^-1 = I
        let i00 = m[0][0] * inv[0][0] + m[0][1] * inv[1][0];
        let i11 = m[1][0] * inv[0][1] + m[1][1] * inv[1][1];
        assert!((i00 - 1.0).abs() < 1e-10);
        assert!((i11 - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_b_matrix() {
        let tri = Triangle::new([
            Point2D::new(0.0, 0.0),
            Point2D::new(1.0, 0.0),
            Point2D::new(0.0, 1.0),
        ]);
        let j = tri.jacobian();
        let j_inv = invert_2x2(&j).unwrap();
        let dn = tri.shape_derivatives(&[0.0, 0.0]);
        let b = compute_b_matrix_2d(&dn, &j_inv);

        assert_eq!(b.len(), 3); // 3 strain components
        assert_eq!(b[0].len(), 6); // 3 nodes * 2 DOFs
    }

    #[test]
    fn test_element_trait() {
        let tri = Triangle::new([
            Point2D::new(0.0, 0.0),
            Point2D::new(1.0, 0.0),
            Point2D::new(0.0, 1.0),
        ]);
        assert_eq!(tri.num_nodes(), 3);
        assert_eq!(tri.dimensions(), 2);

        let tet = Tetrahedron::new([
            Point3D::origin(),
            Point3D::new(1.0, 0.0, 0.0),
            Point3D::new(0.0, 1.0, 0.0),
            Point3D::new(0.0, 0.0, 1.0),
        ]);
        assert_eq!(tet.num_nodes(), 4);
        assert_eq!(tet.dimensions(), 3);
    }
}