#pragma once
/**
 * @file fea_mesh.h
 * @brief 高质量FEA有限元网格生成模块
 *
 * 提供四面体、边界层棱柱、六面体网格生成、自适应细化、
 * 网格质量评估及CAE格式（Abaqus/Ansys/Nastran）导出。
 *
 * ## 功能概览
 * - tetrahedral_mesh: Delaunay 3D + 质量优化四面体网格
 * - boundary_layer_mesh: 沿壁面法向拉伸棱柱边界层
 * - hexahedral_mesh: sweep/extrude 六面体网格
 * - adaptive_refinement: 基于误差估计的局部细化
 * - fea_quality_report: skewness/aspect_ratio/jacobian/orthogonality 报告
 * - export_abaqus / export_ansys / export_nastran: CAE 格式导出
 *
 * @ingroup mesh
 */
#include "vde/mesh/halfedge_mesh.h"
#include "vde/mesh/delaunay_3d.h"
#include "vde/brep/brep.h"
#include <vector>
#include <array>
#include <string>
#include <functional>
#include <cstdio>

namespace vde::mesh {
using core::Point3D;
using core::Vector3D;

// ════════════════════════════════════════════════════════
//  FEA 单元类型
// ════════════════════════════════════════════════════════

/**
 * @brief FEA 单元拓扑类型
 * @ingroup mesh
 */
enum class FEAElementType {
    Tet4,     ///< 4节点线性四面体 (C3D4)
    Tet10,    ///< 10节点二次四面体 (C3D10)
    Hex8,     ///< 8节点线性六面体 (C3D8)
    Hex20,    ///< 20节点二次六面体 (C3D20)
    Wedge6,   ///< 6节点线性三棱柱/楔形 (C3D6)
    Wedge15,  ///< 15节点二次三棱柱 (C3D15)
};

/// 获取单元类型对应的每单元节点数
inline int nodes_per_element(FEAElementType t) {
    switch (t) {
        case FEAElementType::Tet4:   return 4;
        case FEAElementType::Tet10:  return 10;
        case FEAElementType::Hex8:   return 8;
        case FEAElementType::Hex20:  return 20;
        case FEAElementType::Wedge6: return 6;
        case FEAElementType::Wedge15:return 15;
    }
    return 0;
}

/// 获取单元类型名称
inline const char* element_type_name(FEAElementType t) {
    switch (t) {
        case FEAElementType::Tet4:   return "Tet4";
        case FEAElementType::Tet10:  return "Tet10";
        case FEAElementType::Hex8:   return "Hex8";
        case FEAElementType::Hex20:  return "Hex20";
        case FEAElementType::Wedge6: return "Wedge6";
        case FEAElementType::Wedge15:return "Wedge15";
    }
    return "Unknown";
}

// ════════════════════════════════════════════════════════
//  FEAMesh 数据结构
// ════════════════════════════════════════════════════════

/**
 * @brief FEA有限元网格
 *
 * 存储顶点、单元连通性、边界面及单元类型信息。
 * 支持Tet4/Tet10/Hex8/Wedge6等常见FEA单元。
 *
 * @ingroup mesh
 */
struct FEAMesh {
    std::vector<Point3D> vertices;               ///< 顶点坐标数组
    std::vector<std::vector<int>> elements;      ///< 各单元顶点索引（长度 = nodes_per_element）
    FEAElementType element_type = FEAElementType::Tet4;  ///< 单元类型
    std::vector<int> material_ids;               ///< 各单元材料ID（可选，默认0）
    std::vector<std::vector<int>> boundary_faces;///< 边界面（用于施加BC），顶点索引数组
    std::vector<int> boundary_face_tags;         ///< 边界面标签（与boundary_faces一一对应）

    /** @brief 单元总数 */
    [[nodiscard]] size_t num_elements() const { return elements.size(); }
    /** @brief 顶点总数 */
    [[nodiscard]] size_t num_vertices() const { return vertices.size(); }
    /** @brief 边界面总数 */
    [[nodiscard]] size_t num_boundary_faces() const { return boundary_faces.size(); }
    /** @brief 每单元节点数 */
    [[nodiscard]] int npe() const { return nodes_per_element(element_type); }
};

// ════════════════════════════════════════════════════════
//  网格生成参数
// ════════════════════════════════════════════════════════

/**
 * @brief 四面体网格生成参数
 * @ingroup mesh
 */
struct TetMeshParams {
    double min_angle = 15.0;           ///< 最小二面角（度），太小导致刚度矩阵病态
    double max_aspect_ratio = 5.0;     ///< 最大长宽比（外接球半径/最短边）
    double max_size = 0.0;             ///< 最大单元尺寸，0 = 自动（包围盒对角线的 5%）
    double min_size = 0.0;             ///< 最小单元尺寸，0 = max_size/10
    int quality_iterations = 3;        ///< 质量优化迭代次数（Laplacian光顺+边翻转）
    bool keep_surface = true;          ///< 是否保持原表面顶点位置
};

/**
 * @brief 边界层网格参数
 * @ingroup mesh
 */
struct BLPParams {
    double first_cell_height = 0.001;  ///< 第一层棱柱高度（壁面法向距离）
    double growth_rate = 1.2;          ///< 相邻层高度增长率（>1.0）
    int num_layers = 5;                ///< 棱柱层数
    double total_thickness = 0.0;      ///< 边界层总厚度，0 = 自动（first * (growth^n - 1)/(growth - 1)）
    bool smooth_transition = true;     ///< 是否平滑过渡到内部网格
};

/**
 * @brief 六面体网格参数
 * @ingroup mesh
 */
struct HexMeshParams {
    int sweep_layers = 1;              ///< sweep 方向层数
    double sweep_distance = 0.0;       ///< sweep 总距离，0 = 自动（包围盒最小维度）
    bool use_quad_dominant = true;     ///< 源面是否使用四边形主导网格
    double cell_size = 0.0;            ///< 目标单元尺寸，0 = 自动
};

// ════════════════════════════════════════════════════════
//  网格质量评估
// ════════════════════════════════════════════════════════

/**
 * @brief FEA网格质量报告
 *
 * 包含 skewness（偏斜度）、aspect_ratio（长宽比）、
 * jacobian（雅可比）、orthogonality（正交性）四项
 * 核心FEA网格质量指标。
 *
 * @ingroup mesh
 */
struct FEAQualityReport {
    /// 偏斜度 (skewness): 实际单元与理想单元的偏差，[0,1]，0=最优
    double min_skewness = 0.0;
    double max_skewness = 0.0;
    double avg_skewness = 0.0;

    /// 长宽比 (aspect ratio): 最长边/最短边，(Fluent: 1=最佳)
    double min_aspect_ratio = 0.0;
    double max_aspect_ratio = 0.0;
    double avg_aspect_ratio = 0.0;

    /// 归一化雅可比 (scaled jacobian): [0,1]，1=正单元，0=退化，<0=反转
    double min_jacobian = 0.0;
    double max_jacobian = 0.0;
    double avg_jacobian = 0.0;

    /// 正交性 (orthogonality): 面法向与质心连线的夹角余弦，[0,1]，1=完美正交
    double min_orthogonality = 0.0;
    double max_orthogonality = 0.0;
    double avg_orthogonality = 0.0;

    size_t total_elements = 0;         ///< 单元总数
    size_t degenerate_elements = 0;    ///< 退化单元数 (Jacobian < 0.01)
    size_t inverted_elements = 0;      ///< 反转单元数 (Jacobian < 0)

    /// 分级统计 (A=优秀 B=良好 C=合格 D=差 F=不合格，按Jacobian)
    size_t grade_A = 0, grade_B = 0, grade_C = 0, grade_D = 0, grade_F = 0;
};

// ════════════════════════════════════════════════════════
//  自适应细化
// ════════════════════════════════════════════════════════

/**
 * @brief 误差估计函数类型
 *
 * 对给定单元返回估计误差值（≥0），值越大表示该单元越需要细化。
 *
 * @param element_idx 单元索引
 * @param mesh 所属网格
 * @return 误差估计值（非负）
 */
using ErrorEstimator = std::function<double(int element_idx, const FEAMesh& mesh)>;

/**
 * @brief 自适应细化参数
 * @ingroup mesh
 */
struct RefinementParams {
    double error_threshold = 0.5;      ///< 相对误差阈值（高于此值的单元被细化）
    double refinement_fraction = 0.2;  ///< 每次细化单元比例（0~1）
    int max_iterations = 3;            ///< 最大细化迭代次数
    bool smooth_after_refine = true;   ///< 细化后是否光顺
};

// ════════════════════════════════════════════════════════
//  网格生成 API
// ════════════════════════════════════════════════════════

/**
 * @brief 生成四面体FEA网格（Delaunay 3D + 质量优化）
 *
 * 对 BrepModel 进行 tessellation 获取表面三角网格，
 * 采样内部点后进行 Delaunay 3D 四面体化，
 * 再进行 Laplacian 光顺 + 边翻转优化质量。
 *
 * @param body B-Rep实体模型
 * @param params 四面体化参数（最小角、长宽比约束等）
 * @return FEAMesh 四面体网格
 *
 * @code{.cpp}
 * auto box = vde::brep::make_box(1, 1, 1);
 * TetMeshParams p;
 * p.min_angle = 18.0;
 * p.max_aspect_ratio = 4.0;
 * auto mesh = tetrahedral_mesh(box, p);
 * @endcode
 * @ingroup mesh
 */
FEAMesh tetrahedral_mesh(const brep::BrepModel& body,
                         const TetMeshParams& params = {});

/**
 * @brief 生成边界层棱柱网格
 *
 * 沿B-Rep模型壁面法向拉伸多层棱柱单元，
 * 用于CFD边界层捕捉（y+控制）。
 * 先对表面进行三角化，然后每层沿法向偏移。
 *
 * @param body B-Rep实体模型
 * @param params 边界层参数（首层高度、增长率、层数）
 * @return FEAMesh 棱柱网格 (Wedge6)
 *
 * @code{.cpp}
 * BLPParams p;
 * p.first_cell_height = 0.0001;
 * p.growth_rate = 1.15;
 * p.num_layers = 10;
 * auto bl = boundary_layer_mesh(wing, p);
 * @endcode
 * @ingroup mesh
 */
FEAMesh boundary_layer_mesh(const brep::BrepModel& body,
                            const BLPParams& params = {});

/**
 * @brief 生成六面体网格（sweep/extrude策略）
 *
 * 将B-Rep模型的一个面（源面）四边形化，
 * 然后沿法向 sweep 生成六面体层。
 * 适用于拉伸体/棱柱体等 sweepable 几何。
 *
 * @param body B-Rep实体模型
 * @param params 六面体化参数
 * @return FEAMesh 六面体网格 (Hex8)
 *
 * @code{.cpp}
 * auto cylinder = vde::brep::make_cylinder(1, 3);
 * HexMeshParams p;
 * p.sweep_layers = 5;
 * auto mesh = hexahedral_mesh(cylinder, p);
 * @endcode
 * @ingroup mesh
 */
FEAMesh hexahedral_mesh(const brep::BrepModel& body,
                        const HexMeshParams& params = {});

// ════════════════════════════════════════════════════════
//  高级六面体网格生成 API
// ════════════════════════════════════════════════════════

/// 六面体质量优化选项
struct HexOptimizationParams {
    double min_scaled_jacobian = 0.3;    ///< 最小可接受雅可比
    double max_aspect_ratio = 10.0;      ///< 最大长宽比
    int max_iterations = 50;             ///< 最大优化迭代次数
    double convergence_tol = 1e-4;       ///< 收敛容差
    bool untangle = true;                ///< 是否翻转修复（untangle inverted elements）
    bool laplacian_smooth = true;        ///< 是否 Laplacian 光顺
    bool optimization_based_smooth = true; ///< 是否基于优化的光顺
};

/// 六面体质量优化结果
struct HexOptimizationResult {
    FEAMesh optimized_mesh;              ///< 优化后的网格
    FEAQualityReport initial_quality;    ///< 初始质量报告
    FEAQualityReport final_quality;      ///< 最终质量报告
    int iterations = 0;                  ///< 实际迭代次数
    int inverted_fixed = 0;              ///< 修复的反转单元数
    int degenerate_fixed = 0;            ///< 修复的退化单元数
    bool converged = false;              ///< 是否收敛
};

/// 面标识结构
struct FaceRegion {
    std::vector<int> face_indices;       ///< 面的顶点索引列表
    std::string name;                    ///< 面名称
};

/// 多块分区定义
struct BlockRegion {
    std::vector<Point3D> corner_vertices; ///< 8个角顶点（i,j,k 顺序）
    int n_i = 2;                          ///< i方向分段数
    int n_j = 2;                          ///< j方向分段数
    int n_k = 2;                          ///< k方向分段数
};

/**
 * @brief 映射法六面体网格生成
 *
 * 给定 B-Rep 实体的源面和目标面，使用映射法生成结构化六面体网格。
 * 核心思路：将源面四边形网格通过坐标映射变换到目标面，
 * 并在中间层线性插值，形成规则的六面体网格。
 *
 * 算法：
 *  1. 在源面生成四边形网格（重心映射或栅格映射）
 *  2. 在目标面生成对应的四边形网格
 *  3. 沿 sweep 方向线性插值生成中间层
 *  4. 连接各层为六面体单元
 *
 * @param body        B-Rep实体模型
 * @param source_face 源面标识
 * @param target_face 目标面标识
 * @param layers      层数（默认 4）
 * @return FEAMesh 六面体网格 (Hex8)
 *
 * @note 对标 Ansys ICEM CFD Blocking / Pointwise
 * @ingroup mesh
 */
[[nodiscard]] FEAMesh mapped_hex_mesh(
    const brep::BrepModel& body,
    const FaceRegion& source_face,
    const FaceRegion& target_face,
    int layers = 4);

/**
 * @brief 子映射法六面体网格生成
 *
 * 自动将复杂几何体分区为多个可映射子区域，
 * 然后在每个子区域内用映射法生成六面体网格，
 * 最后缝合各子区域网格。
 *
 * 算法：
 *  1. 分析几何体的拓扑特征（角点、棱边）
 *  2. 自动分区为可扫掠（sweepable）子区域
 *  3. 对每个子区域执行映射法六面体化
 *  4. 合并子区域网格，保证界面节点一致
 *
 * @param body B-Rep实体模型
 * @return FEAMesh 六面体网格 (Hex8)
 *
 * @note 对标 Ansys Workbench MultiZone / CUBIT Submap
 * @ingroup mesh
 */
[[nodiscard]] FEAMesh submapped_hex_mesh(
    const brep::BrepModel& body);

/**
 * @brief 多块结构化六面体网格生成
 *
 * 对几何体进行显式多块分解，每个块独立生成结构化六面体网格，
 * 确保块间界面节点一一对应。
 *
 * 算法：
 *  1. 根据用户定义或自动识别的块分区
 *  2. 每个块内使用 transfinite interpolation (TFI) 生成节点
 *  3. 块间共享面节点，保证 conformal 连接
 *  4. 组装块为整体网格
 *
 * @param body   B-Rep实体模型
 * @param blocks 块分区定义列表
 * @return FEAMesh 六面体网格 (Hex8)
 *
 * @note 对标 Ansys ICEM CFD Hexa Blocking / TrueGrid
 * @ingroup mesh
 */
[[nodiscard]] FEAMesh multi_block_hex_mesh(
    const brep::BrepModel& body,
    const std::vector<BlockRegion>& blocks);

/**
 * @brief 六面体网格质量优化与翻转修复
 *
 * 对六面体网格进行质量优化，包括反转变形修复（untangling）、
 * Laplacian 光顺和基于优化的光顺。
 *
 * 算法：
 *  - untangling: 检测负雅可比单元，通过解局部优化问题翻转节点
 *  - Laplacian 光顺: 将内部节点移到邻域质心
 *  - 优化光顺: 最小化目标函数（扭曲度+体积变化）
 *
 * @param mesh   输入 FEA 六面体网格
 * @param params 优化参数
 * @return HexOptimizationResult 优化结果（含优化前后质量报告）
 *
 * @note 对标 CUBIT MeshGems / Pointwise T-Rex Quality
 * @ingroup mesh
 */
[[nodiscard]] HexOptimizationResult hex_quality_optimization(
    const FEAMesh& mesh,
    const HexOptimizationParams& params = {});

// ════════════════════════════════════════════════════════
//  自适应细化
/**
 * 基于误差估计函数对高误差区域进行局部细分。
 * 支持边中点分裂细化（适用于Tet4→4×Tet4子单元）。
 *
 * @param mesh 输入FEA网格
 * @param estimator 误差估计函数
 * @param params 细化参数
 * @return FEAMesh 细化后的网格
 *
 * @code{.cpp}
 * auto estimator = [](int ei, const FEAMesh& m) -> double {
 *     // 基于应力梯度的误差估计
 *     return stress_gradient_norm(ei, m);
 * };
 * auto refined = adaptive_refinement(mesh, estimator);
 * @endcode
 * @ingroup mesh
 */
FEAMesh adaptive_refinement(const FEAMesh& mesh,
                            const ErrorEstimator& estimator,
                            const RefinementParams& params = {});

// ════════════════════════════════════════════════════════
//  质量评估 API
// ════════════════════════════════════════════════════════

/**
 * @brief 评估FEA网格质量
 *
 * 计算四项核心质量指标：
 * - skewness: 偏斜度（实际vs理想单元的偏差）
 * - aspect_ratio: 长宽比
 * - jacobian: 归一化雅可比行列式
 * - orthogonality: 正交性（面法向与质心连线夹角）
 *
 * 自动根据 element_type 选择对应算法（Tet/Hex/Wedge）。
 *
 * @param mesh FEA网格
 * @return FEAQualityReport 质量报告
 *
 * @ingroup mesh
 */
FEAQualityReport fea_quality_report(const FEAMesh& mesh);

// ── 单项质量指标 ──

/**
 * @brief 计算单元的偏斜度
 * @param verts 单元顶点数组
 * @param type 单元类型
 * @return skewness [0,1]，0=理想单元
 */
double element_skewness(const std::vector<Point3D>& verts, FEAElementType type);

/**
 * @brief 计算单元的长宽比
 * @param verts 单元顶点数组
 * @param type 单元类型
 * @return aspect_ratio ≥ 1，1=理想
 */
double element_aspect_ratio(const std::vector<Point3D>& verts, FEAElementType type);

/**
 * @brief 计算单元的归一化雅可比
 * @param verts 单元顶点数组
 * @param type 单元类型
 * @return scaled_jacobian [0,1]，1=正四面体/正方体，0=退化
 */
double element_scaled_jacobian(const std::vector<Point3D>& verts, FEAElementType type);

/**
 * @brief 计算单元的正交性
 * @param verts 单元顶点数组
 * @param type 单元类型
 * @return orthogonality [0,1]，1=完全正交（面法向与对面质心连线平行）
 */
double element_orthogonality(const std::vector<Point3D>& verts, FEAElementType type);

// ════════════════════════════════════════════════════════
//  CAE 格式导出
// ════════════════════════════════════════════════════════

/**
 * @brief 导出 Abaqus 输入文件 (.inp)
 *
 * 生成包含 *NODE 和 *ELEMENT 卡的 Abaqus 输入文件。
 * 支持 Tet4 (C3D4), Tet10 (C3D10), Hex8 (C3D8), Wedge6 (C3D6)。
 *
 * @param mesh FEA网格
 * @param filepath 输出文件路径
 * @return 成功返回 true
 *
 * @ingroup mesh
 */
bool export_abaqus(const FEAMesh& mesh, const std::string& filepath);

/**
 * @brief 导出 Ansys CDB 文件 (.cdb)
 *
 * 生成包含 NBLOCK 和 EBLOCK 的 Ansys CDB 格式文件。
 *
 * @param mesh FEA网格
 * @param filepath 输出文件路径
 * @return 成功返回 true
 *
 * @ingroup mesh
 */
bool export_ansys(const FEAMesh& mesh, const std::string& filepath);

/**
 * @brief 导出 Nastran 批量数据文件 (.bdf/.dat)
 *
 * 生成包含 GRID 和 CTETRA/CHEXA/CPENTA 卡的 Nastran 输入文件。
 *
 * @param mesh FEA网格
 * @param filepath 输出文件路径
 * @return 成功返回 true
 *
 * @ingroup mesh
 */
bool export_nastran(const FEAMesh& mesh, const std::string& filepath);

}  // namespace vde::mesh