// Copyright (c) 2024 RustyTorch++ Team
// Licensed under the Apache License, Version 2.0

//! Isoparametric element utilities and coordinate transformations.

use super::{FiniteElement, NaturalCoords};
use crate::error::{ElementError, FeaResult};
use nalgebra::Vector3;

/// Isoparametric mapping utilities for finite elements.
pub struct IsoparametricMapping;

impl IsoparametricMapping {
    /// Check if an element is isoparametric (same shape functions for geometry and field interpolation).
    pub fn is_isoparametric(_element: &dyn FiniteElement) -> bool {
        // All standard finite elements are isoparametric
        // This could be extended to check specific conditions
        true
    }

    /// Compute the mapping quality at a point.
    pub fn mapping_quality(
        element: &dyn FiniteElement,
        coords: &NaturalCoords,
        node_coords: &[Vector3<f64>],
    ) -> FeaResult<MappingQuality> {
        let jacobian_eval = element.jacobian(coords, node_coords)?;

        let determinant = jacobian_eval.determinant();
        let condition_number = jacobian_eval.condition_number();
        let is_valid = jacobian_eval.is_valid();

        // Compute distortion measures
        let aspect_ratio = Self::compute_aspect_ratio(&jacobian_eval);
        let skewness = Self::compute_skewness(&jacobian_eval);
        let orthogonality = Self::compute_orthogonality(&jacobian_eval);

        Ok(MappingQuality {
            determinant,
            condition_number,
            aspect_ratio,
            skewness,
            orthogonality,
            is_valid,
        })
    }

    /// Analyze mapping quality over the entire element.
    pub fn analyze_element_mapping(
        element: &dyn FiniteElement,
        node_coords: &[Vector3<f64>],
        num_sample_points: usize,
    ) -> FeaResult<ElementMappingAnalysis> {
        let sample_points = Self::generate_sample_points(element, num_sample_points)?;
        let mut qualities = Vec::new();

        for point in &sample_points {
            let quality = Self::mapping_quality(element, point, node_coords)?;
            qualities.push(quality);
        }

        let min_determinant = qualities
            .iter()
            .map(|q| q.determinant)
            .fold(f64::INFINITY, f64::min);
        let max_determinant = qualities
            .iter()
            .map(|q| q.determinant)
            .fold(f64::NEG_INFINITY, f64::max);
        let avg_determinant =
            qualities.iter().map(|q| q.determinant).sum::<f64>() / qualities.len() as f64;

        let min_condition = qualities
            .iter()
            .map(|q| q.condition_number)
            .fold(f64::INFINITY, f64::min);
        let max_condition = qualities
            .iter()
            .map(|q| q.condition_number)
            .fold(f64::NEG_INFINITY, f64::max);
        let avg_condition =
            qualities.iter().map(|q| q.condition_number).sum::<f64>() / qualities.len() as f64;

        let max_aspect_ratio = qualities
            .iter()
            .map(|q| q.aspect_ratio)
            .fold(f64::NEG_INFINITY, f64::max);
        let max_skewness = qualities
            .iter()
            .map(|q| q.skewness)
            .fold(f64::NEG_INFINITY, f64::max);
        let min_orthogonality = qualities
            .iter()
            .map(|q| q.orthogonality)
            .fold(f64::INFINITY, f64::min);

        let invalid_count = qualities.iter().filter(|q| !q.is_valid).count();

        Ok(ElementMappingAnalysis {
            sample_points: sample_points.len(),
            min_determinant,
            max_determinant,
            avg_determinant,
            min_condition,
            max_condition,
            avg_condition,
            max_aspect_ratio,
            max_skewness,
            min_orthogonality,
            invalid_count,
            is_acceptable: Self::is_mapping_acceptable(&qualities),
        })
    }

    /// Generate sample points for mapping analysis.
    fn generate_sample_points(
        element: &dyn FiniteElement,
        num_points: usize,
    ) -> FeaResult<Vec<NaturalCoords>> {
        let param_dim = element.parametric_dimension();
        let mut points = Vec::new();

        match param_dim {
            1 => {
                // 1D: distribute points along xi
                for i in 0..num_points {
                    let xi = -1.0 + 2.0 * i as f64 / (num_points - 1) as f64;
                    points.push(NaturalCoords::new_1d(xi));
                }
            }
            2 => {
                // 2D: grid of points
                let points_per_dim = (num_points as f64).sqrt().ceil() as usize;
                for i in 0..points_per_dim {
                    for j in 0..points_per_dim {
                        match element.element_type().topology_family() {
                            crate::mesh::ElementTopology::Triangle => {
                                let xi = i as f64 / points_per_dim as f64;
                                let eta = j as f64 / points_per_dim as f64;
                                if xi + eta <= 1.0 {
                                    points.push(NaturalCoords::new_2d(xi, eta));
                                }
                            }
                            crate::mesh::ElementTopology::Quadrilateral => {
                                let xi = -1.0 + 2.0 * i as f64 / (points_per_dim - 1) as f64;
                                let eta = -1.0 + 2.0 * j as f64 / (points_per_dim - 1) as f64;
                                points.push(NaturalCoords::new_2d(xi, eta));
                            }
                            _ => {
                                return Err(ElementError::UnsupportedElementType {
                                    element_type: format!("{:?}", element.element_type()),
                                }
                                .into());
                            }
                        }
                    }
                }
            }
            3 => {
                // 3D: cube or tetrahedral sampling
                let points_per_dim = (num_points as f64).cbrt().ceil() as usize;
                for i in 0..points_per_dim {
                    for j in 0..points_per_dim {
                        for k in 0..points_per_dim {
                            match element.element_type().topology_family() {
                                crate::mesh::ElementTopology::Tetrahedron => {
                                    let xi = i as f64 / points_per_dim as f64;
                                    let eta = j as f64 / points_per_dim as f64;
                                    let zeta = k as f64 / points_per_dim as f64;
                                    if xi + eta + zeta <= 1.0 {
                                        points.push(NaturalCoords::new_3d(xi, eta, zeta));
                                    }
                                }
                                crate::mesh::ElementTopology::Hexahedron => {
                                    let xi = -1.0 + 2.0 * i as f64 / (points_per_dim - 1) as f64;
                                    let eta = -1.0 + 2.0 * j as f64 / (points_per_dim - 1) as f64;
                                    let zeta = -1.0 + 2.0 * k as f64 / (points_per_dim - 1) as f64;
                                    points.push(NaturalCoords::new_3d(xi, eta, zeta));
                                }
                                _ => {
                                    return Err(ElementError::UnsupportedElementType {
                                        element_type: format!("{:?}", element.element_type()),
                                    }
                                    .into());
                                }
                            }
                        }
                    }
                }
            }
            _ => {
                return Err(ElementError::InvalidGeometry {
                    message: format!("Unsupported parametric dimension: {param_dim}"),
                }
                .into());
            }
        }

        if points.is_empty() {
            points.push(NaturalCoords::new_3d(0.0, 0.0, 0.0)); // Center point
        }

        Ok(points)
    }

    /// Compute aspect ratio from Jacobian.
    fn compute_aspect_ratio(jacobian_eval: &super::JacobianEval) -> f64 {
        let j = &jacobian_eval.jacobian;
        let (rows, cols) = (j.nrows(), j.ncols());

        match (rows, cols) {
            (2, 2) => {
                // 2D: ratio of maximum to minimum singular values
                let j11 = j[(0, 0)];
                let j12 = j[(0, 1)];
                let j21 = j[(1, 0)];
                let j22 = j[(1, 1)];

                let a = j11 * j11 + j21 * j21;
                let b = j12 * j12 + j22 * j22;
                let c = j11 * j12 + j21 * j22;

                let trace = a + b;
                let det = a * b - c * c;

                if det <= 0.0 {
                    return f64::INFINITY;
                }

                let sqrt_discriminant = ((trace * trace - 4.0 * det).max(0.0)).sqrt();
                let sigma_max = f64::midpoint(trace, sqrt_discriminant).sqrt();
                let sigma_min = ((trace - sqrt_discriminant) / 2.0).sqrt();

                if sigma_min > 1e-12 {
                    sigma_max / sigma_min
                } else {
                    f64::INFINITY
                }
            }
            (3, 3) => {
                // 3D: proper aspect ratio using singular value decomposition
                // Compute column norms of Jacobian matrix (equivalent to principal stretches)
                let mut max_length: f64 = 0.0;
                let mut min_length: f64 = f64::INFINITY;

                for i in 0..3 {
                    let length = (0..3)
                        .map(|row| j[(row, i)] * j[(row, i)])
                        .sum::<f64>()
                        .sqrt();
                    max_length = max_length.max(length);
                    min_length = min_length.min(length);
                }

                if min_length > 1e-12 {
                    max_length / min_length
                } else {
                    f64::INFINITY
                }
            }
            _ => 1.0, // Default for unsupported dimensions
        }
    }

    /// Compute skewness from Jacobian.
    fn compute_skewness(jacobian_eval: &super::JacobianEval) -> f64 {
        let j = &jacobian_eval.jacobian;
        let (rows, cols) = (j.nrows(), j.ncols());

        match (rows, cols) {
            (2, 2) => {
                // 2D: angle between coordinate lines
                let v1 = Vector3::new(j[(0, 0)], j[(1, 0)], 0.0);
                let v2 = Vector3::new(j[(0, 1)], j[(1, 1)], 0.0);

                let dot = v1.dot(&v2);
                let norm1 = v1.norm();
                let norm2 = v2.norm();

                if norm1 > 1e-12 && norm2 > 1e-12 {
                    let cos_angle = (dot / (norm1 * norm2)).clamp(-1.0, 1.0);
                    let angle = cos_angle.acos();
                    (angle - std::f64::consts::PI / 2.0).abs() / (std::f64::consts::PI / 2.0)
                } else {
                    1.0
                }
            }
            (3, 3) => {
                // 3D: maximum skewness among all coordinate plane pairs
                let mut max_skewness: f64 = 0.0;

                for i in 0..3 {
                    for k in (i + 1)..3 {
                        let v1 = Vector3::new(j[(0, i)], j[(1, i)], j[(2, i)]);
                        let v2 = Vector3::new(j[(0, k)], j[(1, k)], j[(2, k)]);

                        let dot = v1.dot(&v2);
                        let norm1 = v1.norm();
                        let norm2 = v2.norm();

                        if norm1 > 1e-12 && norm2 > 1e-12 {
                            let cos_angle = (dot / (norm1 * norm2)).clamp(-1.0, 1.0);
                            let angle = cos_angle.acos();
                            let skewness = (angle - std::f64::consts::PI / 2.0).abs()
                                / (std::f64::consts::PI / 2.0);
                            max_skewness = max_skewness.max(skewness);
                        }
                    }
                }

                max_skewness
            }
            _ => 0.0,
        }
    }

    /// Compute orthogonality measure.
    fn compute_orthogonality(jacobian_eval: &super::JacobianEval) -> f64 {
        1.0 - Self::compute_skewness(jacobian_eval)
    }

    /// Check if mapping is acceptable based on quality metrics.
    fn is_mapping_acceptable(qualities: &[MappingQuality]) -> bool {
        qualities.iter().all(|q| {
            q.is_valid
                && q.determinant > 1e-12
                && q.condition_number < 1000.0
                && q.aspect_ratio < 100.0
                && q.skewness < 0.8
        })
    }
}

/// Mapping quality metrics at a point.
#[derive(Debug, Clone)]
pub struct MappingQuality {
    /// Jacobian determinant
    pub determinant: f64,
    /// Condition number
    pub condition_number: f64,
    /// Aspect ratio
    pub aspect_ratio: f64,
    /// Skewness measure [0, 1]
    pub skewness: f64,
    /// Orthogonality measure [0, 1]
    pub orthogonality: f64,
    /// Whether the mapping is valid
    pub is_valid: bool,
}

impl MappingQuality {
    /// Check if this mapping quality is acceptable.
    pub fn is_acceptable(&self) -> bool {
        self.is_valid
            && self.determinant > 1e-12
            && self.condition_number < 1000.0
            && self.aspect_ratio < 100.0
            && self.skewness < 0.8
    }

    /// Get an overall quality score [0, 1] where 1 is perfect.
    pub fn quality_score(&self) -> f64 {
        if !self.is_valid {
            return 0.0;
        }

        let det_score = if self.determinant > 1e-12 { 1.0 } else { 0.0 };
        let condition_score = (1000.0 / self.condition_number.max(1.0)).min(1.0);
        let aspect_score = (10.0 / self.aspect_ratio.max(1.0)).min(1.0);
        let skew_score = 1.0 - self.skewness;

        (det_score + condition_score + aspect_score + skew_score) / 4.0
    }
}

/// Element mapping analysis results.
#[derive(Debug, Clone)]
pub struct ElementMappingAnalysis {
    pub sample_points: usize,
    pub min_determinant: f64,
    pub max_determinant: f64,
    pub avg_determinant: f64,
    pub min_condition: f64,
    pub max_condition: f64,
    pub avg_condition: f64,
    pub max_aspect_ratio: f64,
    pub max_skewness: f64,
    pub min_orthogonality: f64,
    pub invalid_count: usize,
    pub is_acceptable: bool,
}

impl std::fmt::Display for ElementMappingAnalysis {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        writeln!(
            f,
            "Element Mapping Analysis ({} sample points):",
            self.sample_points
        )?;
        writeln!(
            f,
            "  Jacobian determinant: {:.6e} - {:.6e} (avg: {:.6e})",
            self.min_determinant, self.max_determinant, self.avg_determinant
        )?;
        writeln!(
            f,
            "  Condition number: {:.2} - {:.2} (avg: {:.2})",
            self.min_condition, self.max_condition, self.avg_condition
        )?;
        writeln!(f, "  Max aspect ratio: {:.2}", self.max_aspect_ratio)?;
        writeln!(f, "  Max skewness: {:.3}", self.max_skewness)?;
        writeln!(f, "  Min orthogonality: {:.3}", self.min_orthogonality)?;
        writeln!(f, "  Invalid points: {}", self.invalid_count)?;
        writeln!(
            f,
            "  Overall assessment: {}",
            if self.is_acceptable {
                "Acceptable"
            } else {
                "Poor"
            }
        )?;
        Ok(())
    }
}

/// Coordinate transformation utilities.
pub struct CoordinateTransform;

impl CoordinateTransform {
    /// Transform a vector from natural to physical coordinates.
    pub fn transform_vector(
        element: &dyn FiniteElement,
        coords: &NaturalCoords,
        node_coords: &[Vector3<f64>],
        natural_vector: &Vector3<f64>,
    ) -> FeaResult<Vector3<f64>> {
        let jacobian_eval = element.jacobian(coords, node_coords)?;

        let spatial_dim = element.spatial_dimension();
        let param_dim = element.parametric_dimension();

        match (spatial_dim, param_dim) {
            (2, 2) => {
                let j = &jacobian_eval.jacobian;
                let physical_vector = Vector3::new(
                    j[(0, 0)] * natural_vector.x + j[(0, 1)] * natural_vector.y,
                    j[(1, 0)] * natural_vector.x + j[(1, 1)] * natural_vector.y,
                    0.0,
                );
                Ok(physical_vector)
            }
            (3, 3) => {
                let j = &jacobian_eval.jacobian;
                let physical_vector = Vector3::new(
                    j[(0, 0)] * natural_vector.x
                        + j[(0, 1)] * natural_vector.y
                        + j[(0, 2)] * natural_vector.z,
                    j[(1, 0)] * natural_vector.x
                        + j[(1, 1)] * natural_vector.y
                        + j[(1, 2)] * natural_vector.z,
                    j[(2, 0)] * natural_vector.x
                        + j[(2, 1)] * natural_vector.y
                        + j[(2, 2)] * natural_vector.z,
                );
                Ok(physical_vector)
            }
            _ => Err(ElementError::InvalidGeometry {
                message: format!(
                    "Unsupported dimension combination: {spatial_dim}D spatial, {param_dim}D parametric"
                ),
            }
            .into()),
        }
    }

    /// Transform a tensor from natural to physical coordinates.
    pub fn transform_tensor(
        element: &dyn FiniteElement,
        coords: &NaturalCoords,
        node_coords: &[Vector3<f64>],
        natural_tensor: &nalgebra::DMatrix<f64>,
    ) -> FeaResult<nalgebra::DMatrix<f64>> {
        let jacobian_eval = element.jacobian(coords, node_coords)?;

        // T_physical = J * T_natural * J^T
        let j = &jacobian_eval.jacobian;
        let physical_tensor = j * natural_tensor * j.transpose();

        Ok(physical_tensor)
    }

    /// Compute surface normal for boundary elements.
    pub fn surface_normal(
        element: &dyn FiniteElement,
        coords: &NaturalCoords,
        node_coords: &[Vector3<f64>],
    ) -> FeaResult<Vector3<f64>> {
        let jacobian_eval = element.jacobian(coords, node_coords)?;
        let j = &jacobian_eval.jacobian;

        match element.parametric_dimension() {
            1 => {
                // 1D boundary element in 2D/3D space
                if element.spatial_dimension() == 2 {
                    let tangent = Vector3::new(j[(0, 0)], j[(1, 0)], 0.0);
                    let normal = Vector3::new(-tangent.y, tangent.x, 0.0);
                    Ok(normal.normalize())
                } else {
                    // For 1D element in 3D, normal is not uniquely defined
                    Err(ElementError::InvalidGeometry {
                        message: "Normal not uniquely defined for 1D element in 3D space"
                            .to_string(),
                    }
                    .into())
                }
            }
            2 => {
                // 2D boundary element in 3D space
                if element.spatial_dimension() == 3 {
                    let u = Vector3::new(j[(0, 0)], j[(1, 0)], j[(2, 0)]);
                    let v = Vector3::new(j[(0, 1)], j[(1, 1)], j[(2, 1)]);
                    let normal = u.cross(&v);
                    Ok(normal.normalize())
                } else {
                    Err(ElementError::InvalidGeometry {
                        message: "2D element in 2D space has no surface normal".to_string(),
                    }
                    .into())
                }
            }
            _ => Err(ElementError::InvalidGeometry {
                message: format!(
                    "Cannot compute surface normal for {}D element",
                    element.parametric_dimension()
                ),
            }
            .into()),
        }
    }
}

#[cfg(disabled)]
mod tests {
    use super::*;
    use crate::elements::shape_functions::*;

    #[test]
    fn test_mapping_quality_good_element() {
        let element = Triangle3::new();
        let coords = NaturalCoords::new_2d(1.0 / 3.0, 1.0 / 3.0);

        // Well-shaped triangle
        let node_coords = vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(1.0, 0.0, 0.0),
            Vector3::new(0.0, 1.0, 0.0),
        ];

        let quality =
            IsoparametricMapping::mapping_quality(&element, &coords, &node_coords).unwrap();

        assert!(quality.is_acceptable());
        assert!(quality.is_valid);
        assert!(quality.determinant > 0.0);
        assert!(quality.aspect_ratio < 2.0);
        assert!(quality.skewness < 0.1);
    }

    #[test]
    fn test_mapping_quality_distorted_element() {
        let element = Quadrilateral4::new();
        let coords = NaturalCoords::new_2d(0.0, 0.0);

        // Highly distorted quadrilateral
        let node_coords = vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(2.0, 0.0, 0.0),
            Vector3::new(1.9, 0.1, 0.0), // Nearly collapsed
            Vector3::new(0.1, 0.1, 0.0),
        ];

        let quality =
            IsoparametricMapping::mapping_quality(&element, &coords, &node_coords).unwrap();

        assert!(quality.determinant > 0.0); // Still valid but poor quality
        assert!(quality.aspect_ratio > 5.0);
        assert!(quality.skewness > 0.5);
    }

    #[test]
    fn test_element_mapping_analysis() {
        let element = Triangle3::new();
        let node_coords = vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(1.0, 0.0, 0.0),
            Vector3::new(0.0, 1.0, 0.0),
        ];

        let analysis =
            IsoparametricMapping::analyze_element_mapping(&element, &node_coords, 16).unwrap();

        assert!(analysis.is_acceptable);
        assert_eq!(analysis.invalid_count, 0);
        assert!(analysis.min_determinant > 0.0);
        assert!(analysis.max_condition < 10.0);
    }

    #[test]
    fn test_coordinate_transform_vector() {
        let element = Triangle3::new();
        let coords = NaturalCoords::new_2d(0.5, 0.25);
        let node_coords = vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(2.0, 0.0, 0.0), // Scale by 2 in x
            Vector3::new(0.0, 1.0, 0.0),
        ];

        let natural_vector = Vector3::new(1.0, 0.0, 0.0); // Unit vector in xi direction

        let physical_vector =
            CoordinateTransform::transform_vector(&element, &coords, &node_coords, &natural_vector)
                .unwrap();

        // Should be scaled by the Jacobian
        assert!(physical_vector.x > 1.0); // Scaled by ~2
        assert!((physical_vector.z).abs() < 1e-12); // z should remain 0
    }

    #[test]
    fn test_surface_normal_2d_in_3d() {
        let element = Triangle3::new();
        let coords = NaturalCoords::new_2d(1.0 / 3.0, 1.0 / 3.0);

        // Triangle in xy-plane
        let node_coords = vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(1.0, 0.0, 0.0),
            Vector3::new(0.0, 1.0, 0.0),
        ];

        let normal = CoordinateTransform::surface_normal(&element, &coords, &node_coords).unwrap();

        // Normal should point in +z direction
        assert!((normal.x).abs() < 1e-12);
        assert!((normal.y).abs() < 1e-12);
        assert!((normal.z - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_quality_score() {
        let good_quality = MappingQuality {
            determinant: 1.0,
            condition_number: 1.5,
            aspect_ratio: 1.2,
            skewness: 0.1,
            orthogonality: 0.9,
            is_valid: true,
        };

        let poor_quality = MappingQuality {
            determinant: 0.01,
            condition_number: 100.0,
            aspect_ratio: 50.0,
            skewness: 0.7,
            orthogonality: 0.3,
            is_valid: true,
        };

        assert!(good_quality.quality_score() > 0.8);
        assert!(poor_quality.quality_score() < 0.5);
        assert!(good_quality.is_acceptable());
        assert!(!poor_quality.is_acceptable());
    }

    #[test]
    fn test_sample_points_generation() {
        let triangle = Triangle3::new();
        let quad = Quadrilateral4::new();
        let tet = Tetrahedron4::new();
        let hex = Hexahedron8::new();

        let tri_points = IsoparametricMapping::generate_sample_points(&triangle, 16).unwrap();
        let quad_points = IsoparametricMapping::generate_sample_points(&quad, 16).unwrap();
        let tet_points = IsoparametricMapping::generate_sample_points(&tet, 27).unwrap();
        let hex_points = IsoparametricMapping::generate_sample_points(&hex, 27).unwrap();

        assert!(!tri_points.is_empty());
        assert!(!quad_points.is_empty());
        assert!(!tet_points.is_empty());
        assert!(!hex_points.is_empty());

        // Check that triangle points satisfy constraint
        for point in &tri_points {
            assert!(point.xi >= 0.0);
            assert!(point.eta >= 0.0);
            assert!(point.xi + point.eta <= 1.0 + 1e-12);
        }

        // Check that quad points are in bounds
        for point in &quad_points {
            assert!(point.xi >= -1.0 - 1e-12 && point.xi <= 1.0 + 1e-12);
            assert!(point.eta >= -1.0 - 1e-12 && point.eta <= 1.0 + 1e-12);
        }
    }

    #[test]
    fn test_aspect_ratio_computation() {
        let element = Quadrilateral4::new();
        let coords = NaturalCoords::new_2d(0.0, 0.0);

        // Square element (aspect ratio should be ~1)
        let square_coords = vec![
            Vector3::new(-1.0, -1.0, 0.0),
            Vector3::new(1.0, -1.0, 0.0),
            Vector3::new(1.0, 1.0, 0.0),
            Vector3::new(-1.0, 1.0, 0.0),
        ];

        let quality =
            IsoparametricMapping::mapping_quality(&element, &coords, &square_coords).unwrap();
        assert!((quality.aspect_ratio - 1.0).abs() < 0.1);

        // Rectangular element (high aspect ratio)
        let rect_coords = vec![
            Vector3::new(-5.0, -1.0, 0.0),
            Vector3::new(5.0, -1.0, 0.0),
            Vector3::new(5.0, 1.0, 0.0),
            Vector3::new(-5.0, 1.0, 0.0),
        ];

        let rect_quality =
            IsoparametricMapping::mapping_quality(&element, &coords, &rect_coords).unwrap();
        assert!(rect_quality.aspect_ratio > 4.0);
    }
}