/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "mesh_io.h"

gctl::meshio_element::meshio_element()
{
    enabled = false;
    id = DEFAULT_INVALID_TAG;
    type = NotSet;
}

gctl::meshio_data::meshio_data()
{
    enabled = false;
    d_type = NodeData;
}

void gctl::meshio_data::clear()
{
    enabled = false;
    str_tag.clear();
    real_tag.clear();
    int_tag.clear();
    vert_ptrs.clear();
    elem_ptrs.clear();
    val.clear();
    return;
}

bool gctl::meshio_data::pass_check()
{
    // 检查是否同时连接了顶点和单元体
    if (vert_ptrs.empty() && elem_ptrs.empty()) return false;
    if (!vert_ptrs.empty() && !elem_ptrs.empty()) return false;
    if (str_tag.empty() || real_tag.empty() || int_tag.size() < 3) return false;
    if (!vert_ptrs.empty() && (vert_ptrs.size() != val.size() || int_tag[2] != val.size())) return false;
    if (!elem_ptrs.empty() && (elem_ptrs.size() != val.size() || int_tag[2] != val.size())) return false;
    return true;
}

gctl::meshio_element_group::meshio_element_group()
{
    enabled = false;
    type = NotSet;
    name = "Untitled";
    phys_group = geom_group = part_group = DEFAULT_INVALID_TAG;
}

void gctl::meshio_element_group::enable_elements()
{
    for (size_t e = 0; e < elem_ptrs.size(); e++)
    {
        elem_ptrs[e]->enabled = true;
    }
    return;
}

void gctl::meshio_element_group::disable_elements()
{
    for (size_t e = 0; e < elem_ptrs.size(); e++)
    {
        elem_ptrs[e]->enabled = false;
    }
    return;
}

gctl::mesh_io::mesh_io()
{
    valid_node_size_ = 0;
    valid_elem_size_ = 0;
    valid_group_size_ = 0;
    initialized_ = false;

    elem_gmshcode2type_[0] = NotSet;
    elem_gmshcode2type_[1] = _2NodeLine;
    elem_gmshcode2type_[2] = _3NodeTriangle;
    elem_gmshcode2type_[3] = _4NodeQuadrangle;
    elem_gmshcode2type_[4] = _4NodeTetrahedron;
    elem_gmshcode2type_[5] = _8NodeHexahedron;
    elem_gmshcode2type_[6] = _6NodePrism;
    elem_gmshcode2type_[7] = _5NodePyramid;
    elem_gmshcode2type_[8] = _3NodeSecondOrderLine;
    elem_gmshcode2type_[9] = _6NdoeSecondOrderLine;
    elem_gmshcode2type_[10] = _9NodeSecondOrderQuadrangle;
    elem_gmshcode2type_[11] = _10NodeSecondOrderTetrahedron;
    elem_gmshcode2type_[12] = _27NodeSecondOrderHexahedron;
    elem_gmshcode2type_[13] = _18NodeSecondOrderPrism;
    elem_gmshcode2type_[14] = _14NodeSecondOrderPyramid;
    elem_gmshcode2type_[15] = _1NodePoint;
    elem_gmshcode2type_[16] = _8NodeSecondOrderQuadrangle;
    elem_gmshcode2type_[17] = _20NdoeSecondOrderHexahedron;
    elem_gmshcode2type_[18] = _15NodeSecondOrderPrism;
    elem_gmshcode2type_[19] = _13NodeSecondOrderPyramid;
    elem_gmshcode2type_[20] = _9NodeThirdOrderIncompleteTriangle;
    elem_gmshcode2type_[21] = _10NdoeThirdOrderTriangle;
    elem_gmshcode2type_[22] = _12NodeFourthOrderIncompleteTriangle;
    elem_gmshcode2type_[23] = _15NodeFourthOrderTriangle;
    elem_gmshcode2type_[24] = _15NodeFifthOrderCompleteTriangle;
    elem_gmshcode2type_[25] = _21NodeFifthOrderCompleteTriangle;
    elem_gmshcode2type_[26] = _4NodeThirdOrderEdge;
    elem_gmshcode2type_[27] = _5NodeFourthOrderEdge;
    elem_gmshcode2type_[28] = _6NodeFifthOrderEdge;
    elem_gmshcode2type_[29] = _20NodeThirdOrderTetrahedron;
    elem_gmshcode2type_[30] = _35NodeFourthOrderTetrahedron;
    elem_gmshcode2type_[31] = _56NodeFifithOrderTetrahedron;
    elem_gmshcode2type_[92] = _64NodeThirdOrderHexahedron;
    elem_gmshcode2type_[93] = _125NodeFourthOrderHexahedron;

    elem_vtkcode2type_[3]   = _2NodeLine;
    elem_vtkcode2type_[5]   = _3NodeTriangle;
    elem_vtkcode2type_[9]   = _4NodeQuadrangle;
    elem_vtkcode2type_[10]  = _4NodeTetrahedron;
    elem_vtkcode2type_[12]  = _8NodeHexahedron;
    elem_vtkcode2type_[13]  = _6NodePrism;

    elem_type2gmshcode_[NotSet]                               = 0;
    elem_type2gmshcode_[_2NodeLine]                           = 1;
    elem_type2gmshcode_[_3NodeTriangle]                       = 2;
    elem_type2gmshcode_[_4NodeQuadrangle]                     = 3; 
    elem_type2gmshcode_[_4NodeTetrahedron]                    = 4;
    elem_type2gmshcode_[_8NodeHexahedron]                     = 5;
    elem_type2gmshcode_[_6NodePrism]                          = 6;
    elem_type2gmshcode_[_5NodePyramid]                        = 7;
    elem_type2gmshcode_[_3NodeSecondOrderLine]                = 8;
    elem_type2gmshcode_[_6NdoeSecondOrderLine]                = 9;
    elem_type2gmshcode_[_9NodeSecondOrderQuadrangle]          = 10;
    elem_type2gmshcode_[_10NodeSecondOrderTetrahedron]        = 11;
    elem_type2gmshcode_[_27NodeSecondOrderHexahedron]         = 12;
    elem_type2gmshcode_[_18NodeSecondOrderPrism]              = 13;
    elem_type2gmshcode_[_14NodeSecondOrderPyramid]            = 14;
    elem_type2gmshcode_[_1NodePoint]                          = 15;
    elem_type2gmshcode_[_8NodeSecondOrderQuadrangle]          = 16;
    elem_type2gmshcode_[_20NdoeSecondOrderHexahedron]         = 17;
    elem_type2gmshcode_[_15NodeSecondOrderPrism]              = 18;
    elem_type2gmshcode_[_13NodeSecondOrderPyramid]            = 19;
    elem_type2gmshcode_[_9NodeThirdOrderIncompleteTriangle]   = 20;
    elem_type2gmshcode_[_10NdoeThirdOrderTriangle]            = 21; 
    elem_type2gmshcode_[_12NodeFourthOrderIncompleteTriangle] = 22;
    elem_type2gmshcode_[_15NodeFourthOrderTriangle]           = 23;
    elem_type2gmshcode_[_15NodeFifthOrderCompleteTriangle]    = 24;
    elem_type2gmshcode_[_21NodeFifthOrderCompleteTriangle]    = 25;
    elem_type2gmshcode_[_4NodeThirdOrderEdge]                 = 26;
    elem_type2gmshcode_[_5NodeFourthOrderEdge]                = 27;
    elem_type2gmshcode_[_6NodeFifthOrderEdge]                 = 28;
    elem_type2gmshcode_[_20NodeThirdOrderTetrahedron]         = 29;
    elem_type2gmshcode_[_35NodeFourthOrderTetrahedron]        = 30;
    elem_type2gmshcode_[_56NodeFifithOrderTetrahedron]        = 31;
    elem_type2gmshcode_[_64NodeThirdOrderHexahedron]          = 92;
    elem_type2gmshcode_[_125NodeFourthOrderHexahedron]        = 93;

    elem_type2vtkcode_[_2NodeLine] = 3;
    elem_type2vtkcode_[_3NodeTriangle] = 5;
    elem_type2vtkcode_[_4NodeQuadrangle] = 9;
    elem_type2vtkcode_[_4NodeTetrahedron] = 10;
    elem_type2vtkcode_[_8NodeHexahedron] = 12;
    elem_type2vtkcode_[_6NodePrism] = 13;

    elem_type2size_[_2NodeLine]                           = 2;
    elem_type2size_[_3NodeTriangle]                       = 3;
    elem_type2size_[_4NodeQuadrangle]                     = 4;
    elem_type2size_[_4NodeTetrahedron]                    = 4;
    elem_type2size_[_8NodeHexahedron]                     = 8;
	elem_type2size_[_6NodePrism]                          = 6;
    elem_type2size_[_5NodePyramid]                        = 5; 
    elem_type2size_[_3NodeSecondOrderLine]                = 3; 
    elem_type2size_[_6NdoeSecondOrderLine]                = 6; 
    elem_type2size_[_9NodeSecondOrderQuadrangle]          = 9;
	elem_type2size_[_10NodeSecondOrderTetrahedron]        = 10; 
    elem_type2size_[_27NodeSecondOrderHexahedron]         = 27; 
    elem_type2size_[_18NodeSecondOrderPrism]              = 18; 
    elem_type2size_[_14NodeSecondOrderPyramid]            = 14; 
    elem_type2size_[_1NodePoint]                          = 1;
	elem_type2size_[_8NodeSecondOrderQuadrangle]          = 8; 
    elem_type2size_[_20NdoeSecondOrderHexahedron]         = 20; 
    elem_type2size_[_15NodeSecondOrderPrism]              = 15; 
    elem_type2size_[_13NodeSecondOrderPyramid]            = 13; 
    elem_type2size_[_9NodeThirdOrderIncompleteTriangle]   = 9;
	elem_type2size_[_10NdoeThirdOrderTriangle]            = 10; 
    elem_type2size_[_12NodeFourthOrderIncompleteTriangle] = 12; 
    elem_type2size_[_15NodeFourthOrderTriangle]           = 15; 
    elem_type2size_[_15NodeFifthOrderCompleteTriangle]    = 15; 
    elem_type2size_[_21NodeFifthOrderCompleteTriangle]    = 21;
	elem_type2size_[_4NodeThirdOrderEdge]                 = 4; 
    elem_type2size_[_5NodeFourthOrderEdge]                = 5; 
    elem_type2size_[_6NodeFifthOrderEdge]                 = 6; 
    elem_type2size_[_20NodeThirdOrderTetrahedron]         = 20; 
    elem_type2size_[_35NodeFourthOrderTetrahedron]        = 35;
	elem_type2size_[_56NodeFifithOrderTetrahedron]        = 56; 
    elem_type2size_[_64NodeThirdOrderHexahedron]          = 64; 
    elem_type2size_[_125NodeFourthOrderHexahedron]        = 125;

    elem_type2name_[_2NodeLine]                           = "2-node line";
	elem_type2name_[_3NodeTriangle]                       = "3-node triangle";
	elem_type2name_[_4NodeQuadrangle]                     = "4-node quadrangle";
	elem_type2name_[_4NodeTetrahedron]                    = "4-node tetrahedron";
	elem_type2name_[_8NodeHexahedron]                     = "8-node hexahedron";
	elem_type2name_[_6NodePrism]                          = "6-node prism";
	elem_type2name_[_5NodePyramid]                        = "5-node pyramid";
	elem_type2name_[_3NodeSecondOrderLine]                = "3-node second order line";
	elem_type2name_[_6NdoeSecondOrderLine]                = "6-ndoe second order line";
	elem_type2name_[_9NodeSecondOrderQuadrangle]          = "9-node second order quadrangle";
	elem_type2name_[_10NodeSecondOrderTetrahedron]        = "10-node second order tetrahedron";
	elem_type2name_[_27NodeSecondOrderHexahedron]         = "27-node second order hexahedron";
	elem_type2name_[_18NodeSecondOrderPrism]              = "18-node second order prism";
	elem_type2name_[_14NodeSecondOrderPyramid]            = "14-node second order pyramid";
	elem_type2name_[_1NodePoint]                          = "1-node point";
	elem_type2name_[_8NodeSecondOrderQuadrangle]          = "8-node second order quadrangle";
	elem_type2name_[_20NdoeSecondOrderHexahedron]         = "20-ndoe second order hexahedron";
	elem_type2name_[_15NodeSecondOrderPrism]              = "15-node second order prism";
	elem_type2name_[_13NodeSecondOrderPyramid]            = "13-node second order pyramid";
	elem_type2name_[_9NodeThirdOrderIncompleteTriangle]   = "9-node third order incomplete triangle";
	elem_type2name_[_10NdoeThirdOrderTriangle]            = "10-ndoe third order triangle";
	elem_type2name_[_12NodeFourthOrderIncompleteTriangle] = "12-node fourth order incomplete triangle";
	elem_type2name_[_15NodeFourthOrderTriangle]           = "15-node fourth order triangle";
	elem_type2name_[_15NodeFifthOrderCompleteTriangle]    = "15-node fifth order complete triangle";
	elem_type2name_[_21NodeFifthOrderCompleteTriangle]    = "21-node fifth order complete triangle";
	elem_type2name_[_4NodeThirdOrderEdge]                 = "4-node third order edge";
	elem_type2name_[_5NodeFourthOrderEdge]                = "5-node fourth order edge";
	elem_type2name_[_6NodeFifthOrderEdge]                 = "6-node fifth order edge";
	elem_type2name_[_20NodeThirdOrderTetrahedron]         = "20-node third order tetrahedron";
	elem_type2name_[_35NodeFourthOrderTetrahedron]        = "35-node fourth order tetrahedron";
	elem_type2name_[_56NodeFifithOrderTetrahedron]        = "56-node fifith order tetrahedron";
	elem_type2name_[_64NodeThirdOrderHexahedron]          = "64-node third order hexahedron";
	elem_type2name_[_125NodeFourthOrderHexahedron]        = "125-node fourth order hexahedron";
}
    
gctl::mesh_io::~mesh_io()
{
    reset();
}

void gctl::mesh_io::reset()
{
    if (!nodes_.empty()) nodes_.clear();
    if (!elems_.empty()) elems_.clear();
    if (!nodes_tag_.empty()) nodes_tag_.clear();

    if (!groups_.empty())
    {
        for (size_t i = 0; i < groups_.size(); i++)
        {
            groups_[i].elem_ptrs.clear();
        }
        groups_.clear();
    }

    if (!datas_.empty())
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            datas_[i].clear();
        }
        datas_.clear();
    }

    valid_node_size_ = 0;
    valid_elem_size_ = 0;
    valid_group_size_ = 0;
    initialized_ = false;
}

void gctl::mesh_io::info(std::ostream &ss)
{
    int inttemp,temp_count;
	std::string strtemp;
	std::stringstream stemp;

    ss << "number of nodes: " << valid_node_size_ << ", number of elements: " << valid_elem_size_ << std::endl;
    ss << "name, element type, physical group, geometrical group, partition group, element#" << std::endl;
    for (size_t g = 0; g < groups_.size(); g++)
    {
        if (!groups_[g].enabled) continue;

        ss << groups_[g].name << ", ";
        ss << elem_name(groups_[g].type) << ", ";
        ss << groups_[g].phys_group << ", ";
        ss << groups_[g].geom_group << ", ";
        ss << groups_[g].part_group << ", ";
        ss << groups_[g].elem_ptrs.size() << std::endl;
    }

    size_t d_size;
    double min, max;
    for (size_t d = 0; d < datas_.size(); d++)
    {
        if (!datas_[d].enabled) continue;

        d_size = 0;
        min = GCTL_BDL_MAX;
        max = GCTL_BDL_MIN;

        if (datas_[d].d_type == NodeData)
        {
            ss << "nodedata: \"";
            for (size_t l = 0; l < datas_[d].val.size(); l++)
            {
                if (datas_[d].vert_ptrs[l]->id != DEFAULT_INVALID_TAG)
                {
                    min = std::min(min, datas_[d].val[l]);
                    max = std::max(max, datas_[d].val[l]);
                    d_size++;
                }
            }
        }
        else
        {
            ss << "elementdata: \"";
            for (size_t l = 0; l < datas_[d].val.size(); l++)
            {
                if (datas_[d].elem_ptrs[l]->enabled)
                {
                    min = std::min(min, datas_[d].val[l]);
                    max = std::max(max, datas_[d].val[l]);
                    d_size++;
                }
            }
        }

        if (d_size) ss << datas_[d].str_tag[0] << "\", data number: " << d_size << ", mini: " << min << ", maxi = " << max << std::endl;
        else ss << datas_[d].str_tag[0] << "\", no valid value found.\n";
    }
    
    return;
}

void gctl::mesh_io::edit_data(switch_type_e swt, std::string dataname)
{
    if (dataname == "null" && swt == Enable)
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            datas_[i].enabled = true;
        }
    }
    else if (dataname == "null" && swt == Disable)
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            datas_[i].enabled = false;
        }
    }
    else
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            if (datas_[i].str_tag.front() == dataname && swt == Enable) datas_[i].enabled = true;
            if (datas_[i].str_tag.front() == dataname && swt == Disable) datas_[i].enabled = false;
        }
    }
    return;
}

void gctl::mesh_io::edit_group(switch_type_e swt, element_type_enum e_type)
{
    for (size_t g = 0; g < groups_.size(); g++)
    {
        if ((e_type == NotSet || groups_[g].type == e_type) && swt == Enable)
        {
            groups_[g].enabled = true;
            groups_[g].enable_elements();
        }
        else if ((e_type == NotSet || groups_[g].type == e_type) && swt == Disable)
        {
            groups_[g].enabled = false;
            groups_[g].disable_elements();
        }
    }

    update_indexing();
    return;
}

void gctl::mesh_io::edit_group(switch_type_e swt, std::string grp_name)
{
    for (size_t g = 0; g < groups_.size(); g++)
    {
        if (groups_[g].name == grp_name && swt == Enable)
        {
            groups_[g].enabled = true;
            groups_[g].enable_elements();
        }
        else if (groups_[g].name == grp_name && swt == Disable)
        {
            groups_[g].enabled = false;
            groups_[g].disable_elements();
        }
    }

    update_indexing();
    return;
}

void gctl::mesh_io::edit_group(switch_type_e swt, element_tag_enum tag_type, int tag)
{
    for (size_t g = 0; g < groups_.size(); g++)
    {
        if (((tag_type == PhysicalTag && groups_[g].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[g].geom_group == tag) || 
            (tag_type == PartitionTag && groups_[g].part_group == tag)) && 
            swt == Enable)
        {
            groups_[g].enabled = true;
            groups_[g].enable_elements();
        }
        else if (((tag_type == PhysicalTag && groups_[g].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[g].geom_group == tag) || 
            (tag_type == PartitionTag && groups_[g].part_group == tag)) && 
            swt == Disable)
        {
            groups_[g].enabled = false;
            groups_[g].disable_elements();
        }
    }

    update_indexing();
    return;
}

void gctl::mesh_io::edit_group(element_tag_enum anchor_type, int anchor_group, std::string new_name)
{
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((anchor_type == PhysicalTag && groups_[i].phys_group == anchor_group) ||
            (anchor_type == GeometryTag && groups_[i].geom_group == anchor_group) ||
            (anchor_type == PartitionTag && groups_[i].part_group == anchor_group))
        {
            groups_[i].name = new_name;
        }
    }

    // No need to update indexing here.
    return;
}

void gctl::mesh_io::edit_group(element_tag_enum anchor_type, int anchor_group, element_tag_enum tar_type, int tar_group)
{
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((anchor_type == PhysicalTag && groups_[i].phys_group == anchor_group) ||
            (anchor_type == GeometryTag && groups_[i].geom_group == anchor_group) ||
            (anchor_type == PartitionTag && groups_[i].part_group == anchor_group))
        {
            if (tar_type == PhysicalTag) groups_[i].phys_group = tar_group;
            if (tar_type == GeometryTag) groups_[i].geom_group = tar_group;
            if (tar_type == PartitionTag) groups_[i].part_group = tar_group;
        }
    }

    // 检查是否与现有分组重叠 有则合并
    sort_groups();
    return;
}

const gctl::array<int> &gctl::mesh_io::get_node_tag()
{
    return nodes_tag_;
}

int gctl::mesh_io::get_tag(element_tag_enum anchor_type, std::string anchor_name)
{
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (anchor_type == PhysicalTag && groups_[i].name == anchor_name) return groups_[i].phys_group;
        if (anchor_type == GeometryTag && groups_[i].name == anchor_name) return groups_[i].geom_group;
        if (anchor_type == PartitionTag && groups_[i].name == anchor_name) return groups_[i].part_group;
    }
    return DEFAULT_INVALID_TAG;
}

const gctl::array<gctl::vertex3dc> &gctl::mesh_io::get_nodes()
{
    return nodes_;
}

const gctl::array<gctl::meshio_element> &gctl::mesh_io::get_elems()
{
    return elems_;
}

size_t gctl::mesh_io::element_size(element_type_enum e_type)
{
    size_t s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled)
        {
            s += groups_[i].elem_ptrs.size();
        }
    }
    return s;
}

size_t gctl::mesh_io::element_size(element_tag_enum tag_type, int tag)
{
    size_t s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[i].geom_group == tag) || 
            (tag_type == PartitionTag && groups_[i].part_group == tag) && 
            groups_[i].enabled)
        {
            s += groups_[i].elem_ptrs.size();
        }
    }
    return s;
}

size_t gctl::mesh_io::element_size(std::string phys_name)
{
    size_t s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].name == phys_name)
        {
            s += groups_[i].elem_ptrs.size();
        }
    }
    return s;
}

void gctl::mesh_io::convert_tags_to_data(element_tag_enum tag_type)
{
    meshio_data tmp_data;
    if (tag_type == NodeTag && (!nodes_tag_.empty()))
    {
        tmp_data.enabled = true;
        tmp_data.d_type = NodeData;
        tmp_data.str_tag.resize(1, "Node Tag");
        tmp_data.real_tag.resize(1, 0.0);
        tmp_data.int_tag.resize(3, 0);
        tmp_data.int_tag[1] = 1;
        tmp_data.int_tag[2] = valid_node_size_;

        tmp_data.val.resize(valid_node_size_);
        tmp_data.vert_ptrs.resize(valid_node_size_);

        size_t c = 0;
        for (size_t i = 0; i < nodes_.size(); i++)
        {
            if (nodes_[i].id != DEFAULT_INVALID_TAG)
            {
                tmp_data.val[c] = (double) nodes_tag_[i];
                tmp_data.vert_ptrs[c] = nodes_.get(i);
                c++;
            }
        }

        datas_.push_back(tmp_data);
    }
    else
    {
        tmp_data.enabled = true;
        tmp_data.d_type = ElemData;

        if (tag_type == PhysicalTag) tmp_data.str_tag.resize(1, "Physical Tag");
        else if (tag_type == GeometryTag) tmp_data.str_tag.resize(1, "Geometry Tag");
        else if (tag_type == PartitionTag) tmp_data.str_tag.resize(1, "Partition Tag");

        tmp_data.real_tag.resize(1, 0.0);
        tmp_data.int_tag.resize(3, 0);
        tmp_data.int_tag[1] = 1;
        tmp_data.int_tag[2] = valid_elem_size_;

        tmp_data.val.resize(valid_elem_size_);
        tmp_data.elem_ptrs.resize(valid_elem_size_);

        size_t c = 0;
        for (size_t g = 0; g < groups_.size(); g++)
        {
            if (groups_[g].enabled)
            {
                for (size_t e = 0; e < groups_[g].elem_ptrs.size(); e++)
                {
                    if (tag_type == PhysicalTag) tmp_data.val[c] = (double) groups_[g].phys_group;
                    if (tag_type == GeometryTag) tmp_data.val[c] = (double) groups_[g].geom_group;
                    if (tag_type == PartitionTag) tmp_data.val[c] = (double) groups_[g].part_group;
                    tmp_data.elem_ptrs[c] = groups_[g].elem_ptrs[e];
                    c++;
                }
            }
        }

        datas_.push_back(tmp_data);
    }
    return;
}

void gctl::mesh_io::export_elements_to(array<triangle> &tris, std::string phys_name)
{
    size_t s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && 
            (groups_[i].name == phys_name || phys_name == "All"))
        {
            s += groups_[i].elem_ptrs.size();
        }
    }
    
    tris.resize(s);

    s = 0;
    for (size_t i = 0; i < elems_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && 
            (groups_[i].name == phys_name || phys_name == "All"))
        {
            for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++)
            {
                tris[s].id = groups_[i].elem_ptrs[e]->id;
                tris[s].vert[0] = groups_[i].elem_ptrs[e]->vert_ptrs[0];
                tris[s].vert[1] = groups_[i].elem_ptrs[e]->vert_ptrs[1];
                tris[s].vert[2] = groups_[i].elem_ptrs[e]->vert_ptrs[2];
                s++;
            }
        }
    }
    return;
}

void gctl::mesh_io::export_elements_to(array<triangle> &tris, element_tag_enum tag_type, int tag)
{
    size_t s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && 
            ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[i].geom_group == tag)  || 
            (tag_type == PartitionTag && groups_[i].part_group == tag)))
        {
            s += groups_[i].elem_ptrs.size();
        }
    }
    
    tris.resize(s);

    s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && 
            ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[i].geom_group == tag)  || 
            (tag_type == PartitionTag && groups_[i].part_group == tag)))
        {
            for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++)
            {
                tris[s].id = groups_[i].elem_ptrs[e]->id;
                tris[s].vert[0] = groups_[i].elem_ptrs[e]->vert_ptrs[0];
                tris[s].vert[1] = groups_[i].elem_ptrs[e]->vert_ptrs[1];
                tris[s].vert[2] = groups_[i].elem_ptrs[e]->vert_ptrs[2];
                s++;
            }
        }
    }
    return;
}

void gctl::mesh_io::export_elements_to(array<tetrahedron> &tets, std::string phys_name)
{
    size_t s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && 
            (groups_[i].name == phys_name || phys_name == "All"))
        {
            s += groups_[i].elem_ptrs.size();
        }
    }
    
    tets.resize(s);

    s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && 
            (groups_[i].name == phys_name || phys_name == "All"))
        {
            for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++)
            {
                tets[s].set(groups_[i].elem_ptrs[e]->vert_ptrs[0], 
                    groups_[i].elem_ptrs[e]->vert_ptrs[1], 
                    groups_[i].elem_ptrs[e]->vert_ptrs[2], 
                    groups_[i].elem_ptrs[e]->vert_ptrs[3], 
                    groups_[i].elem_ptrs[e]->id);
                s++;
            }
        }
    }
    return;
}

void gctl::mesh_io::export_elements_to(array<tetrahedron> &tets, element_tag_enum tag_type, int tag)
{
    size_t s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && 
            (tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[i].geom_group == tag) || 
            (tag_type == PartitionTag && groups_[i].part_group == tag))
        {
            s += groups_[i].elem_ptrs.size();
        }
    }
    
    tets.resize(s);

    s = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && 
            (tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[i].geom_group == tag) || 
            (tag_type == PartitionTag && groups_[i].part_group == tag))
        {
            for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++)
            {
                tets[s].set(groups_[i].elem_ptrs[e]->vert_ptrs[0], 
                    groups_[i].elem_ptrs[e]->vert_ptrs[1], 
                    groups_[i].elem_ptrs[e]->vert_ptrs[2], 
                    groups_[i].elem_ptrs[e]->vert_ptrs[3], 
                    groups_[i].elem_ptrs[e]->id);
                s++;
            }
        }
    }
    return;
}

void gctl::mesh_io::get_gmsh_physical_groups(std::vector<gmsh_physical_group> &g_groups)
{
    if (!g_groups.empty()) g_groups.clear();

    gmsh_physical_group tmp_group;
    bool not_found;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (g_groups.empty() && groups_[i].enabled)
        {
            tmp_group.name = groups_[i].name;
            tmp_group.phys_tag = groups_[i].phys_group;
            tmp_group.dim_tag = groups_[i].part_group;
            g_groups.push_back(tmp_group);
        }
        else
        {
            not_found = true;
            for (size_t g = 0; g < g_groups.size(); g++)
            {
                if (groups_[i].part_group == g_groups[g].dim_tag &&
                    groups_[i].phys_group == g_groups[g].phys_tag)
                {
                    not_found = false;
                    break;
                }
            }

            if (not_found && groups_[i].enabled)
            {
                tmp_group.name = groups_[i].name;
                tmp_group.phys_tag = groups_[i].phys_group;
                tmp_group.dim_tag = groups_[i].part_group;
                g_groups.push_back(tmp_group);
            }
        }
    }
    return;
}

int gctl::mesh_io::if_saved_data(std::string name, mesh_data_type_e type)
{
    for (size_t i = 0; i < datas_.size(); i++)
    {
        if (datas_[i].str_tag.front() == name &&
            datas_[i].d_type == type) return i;
    }
    return -1;
}

gctl::meshio_data &gctl::mesh_io::get_data(std::string name, mesh_data_type_e type)
{
    int id = if_saved_data(name, type);
    if (id == -1) throw std::runtime_error("[gctl::mesh_io::get_data] Data not found.");
    return datas_[id];
}

gctl::meshio_data *gctl::mesh_io::get_data_ptr(std::string name, mesh_data_type_e type)
{
    int id = if_saved_data(name, type);
    if (id == -1) throw std::runtime_error("[gctl::mesh_io::get_data] Data not found.");
    return &datas_[id];
}

void gctl::mesh_io::add_node_data(std::string name, const array<double> &data)
{
    size_t s = nodes_.size();
    if (data.size()!= s) throw std::runtime_error("[gctl::mesh_io::create_node_data] Incompatible data size.");

    int d_id = if_saved_data(name, NodeData);
    if (d_id != -1)
    {
        for (size_t i = 0; i < s; i++)
        {
            datas_[d_id].val[i] = data[i];
        }
        return;
    }

    meshio_data new_data;
    new_data.enabled = true;
    new_data.d_type = NodeData;
    new_data.str_tag.resize(1, name);
    new_data.real_tag.resize(1, 0.0);
    new_data.int_tag.resize(3, 0);
    new_data.int_tag[1] = 1;
    new_data.int_tag[2] = s;
    new_data.vert_ptrs.resize(s, nullptr);
    new_data.val.resize(s);

    for (size_t i = 0; i < s; i++)
    {
        new_data.vert_ptrs[i] = nodes_.get(i);
        new_data.val[i] = data[i];
    }

    datas_.push_back(new_data);
    return;
}

void gctl::mesh_io::add_node_data(std::string name, const array<double> &data, const array<bool> &boolen)
{
    size_t s = nodes_.size();
    if (data.size()!= s || boolen.size() != s) throw std::runtime_error("[gctl::mesh_io::create_node_data] Incompatible data size.");

    s = 0;
    for (size_t i = 0; i < nodes_.size(); i++)
    {
        if (boolen[i]) s++;
    }

    int d_id = if_saved_data(name, NodeData);
    if (d_id != -1)
    {
        datas_[d_id].val.resize(s);
        datas_[d_id].vert_ptrs.resize(s, nullptr);
        datas_[d_id].int_tag[2] = s;

        s = 0;
        for (size_t i = 0; i < nodes_.size(); i++)
        {
            if (boolen[i])
            {
                datas_[d_id].vert_ptrs[s] = nodes_.get(i);
                datas_[d_id].val[s] = data[i];
                s++;
            }
        }
        return;
    }

    meshio_data new_data;
    new_data.enabled = true;
    new_data.d_type = NodeData;
    new_data.str_tag.resize(1, name);
    new_data.real_tag.resize(1, 0.0);
    new_data.int_tag.resize(3, 0);
    new_data.int_tag[1] = 1;
    new_data.int_tag[2] = s;
    new_data.vert_ptrs.resize(s, nullptr);
    new_data.val.resize(s);

    s = 0;
    for (size_t i = 0; i < nodes_.size(); i++)
    {
        if (boolen[i])
        {
            new_data.vert_ptrs[s] = nodes_.get(i);
            new_data.val[s] = data[i];
            s++;
        }
    }

    datas_.push_back(new_data);
    return;
}

void gctl::mesh_io::add_element_data(std::string name, const array<double> &data, element_type_enum e_type)
{
    size_t e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled) e += groups_[i].elem_ptrs.size();
    }

    if (data.size() != e) throw std::runtime_error("[gctl::mesh_io::create_element_data] Incompatible data size.");

    int d_id = if_saved_data(name, ElemData);
    if (d_id != -1)
    {
        datas_[d_id].int_tag[2] = e;
        datas_[d_id].elem_ptrs.resize(e, nullptr);
        datas_[d_id].val.resize(e, 0.0);

        e = 0;
        for (size_t i = 0; i < groups_.size(); i++)
        {
            if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled)
            {
                for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
                {
                    datas_[d_id].elem_ptrs[e] = groups_[i].elem_ptrs[j];
                    datas_[d_id].val[e] = data[e];
                    e++;
                }
            }
        }
        return;
    }

    meshio_data new_data;
    new_data.enabled = true;
    new_data.d_type = ElemData;
    new_data.str_tag.resize(1, name);
    new_data.real_tag.resize(1, 0.0);
    new_data.int_tag.resize(3, 0);
    new_data.int_tag[1] = 1;
    new_data.int_tag[2] = e;
    new_data.elem_ptrs.resize(e, nullptr);
    new_data.val.resize(e, 0.0);

    e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled)
        {
            for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
            {
                new_data.elem_ptrs[e] = groups_[i].elem_ptrs[j];
                new_data.val[e] = data[e];
                e++;
            }
        }
    }

    datas_.push_back(new_data);
    return;
}

void gctl::mesh_io::add_element_data(std::string name, const array<double> &data, element_tag_enum tag_type, int tag)
{
    size_t e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[i].geom_group == tag) || 
            (tag_type == PartitionTag && groups_[i].part_group == tag) &&
            groups_[i].enabled)
        {
            e += groups_[i].elem_ptrs.size();
        }
    }

    if (data.size() != e) throw std::runtime_error("[gctl::mesh_io::create_element_data] Incompatible data size.");

    int d_id = if_saved_data(name, ElemData);
    if (d_id != -1) 
    {
        array<meshio_element*> more_elem_ptrs(e, nullptr);
        array<double> more_val(e, 0.0);

        e = 0;
        for (size_t i = 0; i < groups_.size(); i++)
        {
            if ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
                (tag_type == GeometryTag && groups_[i].geom_group == tag) || 
                (tag_type == PartitionTag && groups_[i].part_group == tag) &&
                groups_[i].enabled)
            {
                for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
                {
                    more_elem_ptrs[e] = groups_[i].elem_ptrs[j];
                    more_val[e] = data[e];
                    e++;
                }
            }
        }

        datas_[d_id].int_tag[2] += e;
        datas_[d_id].elem_ptrs.concat(more_elem_ptrs);
        datas_[d_id].val.concat(more_val);
        return;
    }

    meshio_data new_data;
    new_data.enabled = true;
    new_data.d_type = ElemData;
    new_data.str_tag.resize(1, name);
    new_data.real_tag.resize(1, 0.0);
    new_data.int_tag.resize(3, 0);
    new_data.int_tag[1] = 1;
    new_data.int_tag[2] = e;
    new_data.elem_ptrs.resize(e, nullptr);
    new_data.val.resize(e, 0.0);

    e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || 
            (tag_type == GeometryTag && groups_[i].geom_group == tag) || 
            (tag_type == PartitionTag && groups_[i].part_group == tag) && 
            groups_[i].enabled)
        {
            for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
            {
                new_data.elem_ptrs[e] = groups_[i].elem_ptrs[j];
                new_data.val[e] = data[e];
                e++;
            }
        }
    }

    datas_.push_back(new_data);
    return;
}

void gctl::mesh_io::add_element_data(std::string name, std::string phys_name, const array<double> &data)
{
    size_t e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].name == phys_name)
        {
            e += groups_[i].elem_ptrs.size();
        }
    }

    if (data.size() != e) throw std::runtime_error("[gctl::mesh_io::create_element_data] Incompatible data size.");

    int d_id = if_saved_data(name, ElemData);
    if (d_id != -1) 
    {
        array<meshio_element*> more_elem_ptrs(e, nullptr);
        array<double> more_val(e, 0.0);

        e = 0;
        for (size_t i = 0; i < groups_.size(); i++)
        {
            if (groups_[i].enabled && groups_[i].name == phys_name)
            {
                for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
                {
                    more_elem_ptrs[e] = groups_[i].elem_ptrs[j];
                    more_val[e] = data[e];
                    e++;
                }
            }
        }

        datas_[d_id].int_tag[2] += e;
        datas_[d_id].elem_ptrs.concat(more_elem_ptrs);
        datas_[d_id].val.concat(more_val);
        return;
    }
    
    meshio_data new_data;
    new_data.enabled = true;
    new_data.d_type = ElemData;
    new_data.str_tag.resize(1, name);
    new_data.real_tag.resize(1, 0.0);
    new_data.int_tag.resize(3, 0);
    new_data.int_tag[1] = 1;
    new_data.int_tag[2] = e;
    new_data.elem_ptrs.resize(e, nullptr);
    new_data.val.resize(e, 0.0);

    e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].name == phys_name)
        {
            for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
            {
                new_data.elem_ptrs[e] = groups_[i].elem_ptrs[j];
                new_data.val[e] = data[e];
                e++;
            }
        }
    }

    datas_.push_back(new_data);
    return;
}

void gctl::mesh_io::add_element_data(std::string name, std::string phys_name, double phys_val)
{
    size_t e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].name == phys_name)
        {
            e += groups_[i].elem_ptrs.size();
        }
    }

    int d_id = if_saved_data(name, ElemData);
    if (d_id != -1) 
    {
        array<meshio_element*> more_elem_ptrs(e, nullptr);
        array<double> more_val(e, 0.0);

        e = 0;
        for (size_t i = 0; i < groups_.size(); i++)
        {
            if (groups_[i].enabled && groups_[i].name == phys_name)
            {
                for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
                {
                    more_elem_ptrs[e] = groups_[i].elem_ptrs[j];
                    more_val[e] = phys_val;
                    e++;
                }
            }
        }

        datas_[d_id].int_tag[2] += e;
        datas_[d_id].elem_ptrs.concat(more_elem_ptrs);
        datas_[d_id].val.concat(more_val);
        return;
    }
    
    meshio_data new_data;
    new_data.enabled = true;
    new_data.d_type = ElemData;
    new_data.str_tag.resize(1, name);
    new_data.real_tag.resize(1, 0.0);
    new_data.int_tag.resize(3, 0);
    new_data.int_tag[1] = 1;
    new_data.int_tag[2] = e;
    new_data.elem_ptrs.resize(e, nullptr);
    new_data.val.resize(e, 0.0);

    e = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled && groups_[i].name == phys_name)
        {
            for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++)
            {
                new_data.elem_ptrs[e] = groups_[i].elem_ptrs[j];
                new_data.val[e] = phys_val;
                e++;
            }
        }
    }

    datas_.push_back(new_data);
    return;
}

void gctl::mesh_io::read_triangle_ascii(std::string filename, index_packed_e is_packed)
{
    if (initialized_ == true)
    {
        throw std::runtime_error("[gctl::mesh_io::read_triangle_ascii()] The mesh is already initialized. Run clear() first.");
    }

    array<vertex2dc> node_2d;
    array<int> node_tag;
    read_Triangle_node(filename, node_2d, is_packed, &node_tag);
    if (!node_tag.empty()) nodes_tag_ = node_tag;

    array<triangle2d> tri_2d;
    matrix<double> tri_attri;
    read_Triangle_element(filename, tri_2d, node_2d, is_packed, &tri_attri);

    valid_node_size_ = node_2d.size();
    valid_elem_size_ = tri_2d.size();
    valid_group_size_ = 0;

    nodes_.resize(valid_node_size_);
    for (size_t i = 0; i < valid_node_size_; i++)
    {
        nodes_[i].id = i;
        nodes_[i].x = node_2d[i].x;
        nodes_[i].y = node_2d[i].y;
        nodes_[i].z = 0.0;
    }

    elems_.resize(valid_elem_size_);
    matrix<int> int_tag(valid_elem_size_, 3, DEFAULT_INVALID_TAG);

    for (size_t i = 0; i < valid_elem_size_; i++)
    {
        elems_[i].id = i;
        elems_[i].neigh_ptrs.resize(3, nullptr);
        elems_[i].type = _3NodeTriangle;

        int_tag[i][2] = 2;

        elems_[i].vert_ptrs.resize(3);
        elems_[i].vert_ptrs[0] = nodes_.get(tri_2d[i].vert[0]->id);
        elems_[i].vert_ptrs[1] = nodes_.get(tri_2d[i].vert[1]->id);
        elems_[i].vert_ptrs[2] = nodes_.get(tri_2d[i].vert[2]->id);
    }

    if (!tri_attri.empty())
    {
        for (size_t i = 0; i < valid_elem_size_; i++)
        {
            int_tag[i][1] = (int) tri_attri[i][0];
        }
    }
    
    // 整理单元体组
    meshio_element_group tmp_group;
    tmp_group.type = elems_[0].type;
    tmp_group.phys_group = int_tag[0][0];
    tmp_group.geom_group = int_tag[0][1];
    tmp_group.part_group = int_tag[0][2];
    tmp_group.elem_ptrs.push_back(elems_.get(0));
    groups_.push_back(tmp_group);

    bool not_found;
    for (size_t i = 1; i < elems_.size(); i++)
    {
        not_found = true;
        for (size_t g = 0; g < groups_.size(); g++)
        {
            if (groups_[g].type == elems_[i].type &&
                groups_[g].phys_group == int_tag[i][0] &&
                groups_[g].geom_group == int_tag[i][1] &&
                groups_[g].part_group == int_tag[i][2])
            {
                groups_[g].elem_ptrs.push_back(elems_.get(i));
                not_found = false;
                break;
            }
        }
        
        if (not_found)
        {
            tmp_group.elem_ptrs.clear();
            tmp_group.type = elems_[i].type;
            tmp_group.phys_group = int_tag[i][0];
            tmp_group.geom_group = int_tag[i][1];
            tmp_group.part_group = int_tag[i][2];
            tmp_group.elem_ptrs.push_back(elems_.get(i));
            groups_.push_back(tmp_group);
        }
    }

    std::string n_file = filename + ".neigh";
    if (!access(n_file.c_str(), R_OK))
    {
        read_Triangle_neighbor(filename, tri_2d, is_packed);

        for (size_t i = 0; i < valid_elem_size_; i++)
        {
            for (size_t j = 0; j < 3; j++)
            {
                if (tri_2d[i].neigh[j] != nullptr) elems_[i].neigh_ptrs[j] = elems_.get(tri_2d[i].neigh[j]->id);
            }
        }
    }

     // 最后将所有的组别和单元体设置为可用
    valid_group_size_ = groups_.size();
    for (size_t g = 0; g < groups_.size(); g++)
    {
        groups_[g].enabled = true;
        groups_[g].enable_elements();
    }
    initialized_ = true;
    return;
}

void gctl::mesh_io::read_tetgen_ascii(std::string filename, index_packed_e is_packed)
{
    if (initialized_ == true)
    {
        throw std::runtime_error("[gctl::mesh_io::read_tetgen_ascii()] The mesh is already initialized. Run clear() first.");
    }

    array<vertex3dc> node_3d;
    array<int> node_tag;
    read_Tetgen_node(filename, node_3d, is_packed, &node_tag);
    if (!node_tag.empty()) nodes_tag_ = node_tag;

    array<tetrahedron> tet_3d;
    array<int> tet_tag;
    read_Tetgen_element(filename, tet_3d, node_3d, is_packed, &tet_tag);

    std::string f_file = filename + ".face";
    array<triangle> tri_3d;
    array<int> tri_tag;
    if (!access(f_file.c_str(), R_OK))
    {    
        read_Tetgen_face(filename, tri_3d, node_3d, is_packed, &tri_tag);
    }

    valid_node_size_ = node_3d.size();
    valid_elem_size_ = tet_3d.size() + tri_3d.size();
    valid_group_size_ = 0;

    nodes_.resize(valid_node_size_);
    for (size_t i = 0; i < valid_node_size_; i++)
    {
        nodes_[i].id = i;
        nodes_[i].x = node_3d[i].x;
        nodes_[i].y = node_3d[i].y;
        nodes_[i].z = node_3d[i].z;
    }

    elems_.resize(valid_elem_size_);
    matrix<int> int_tag(valid_elem_size_, 3, DEFAULT_INVALID_TAG);

    for (size_t i = 0; i < tet_3d.size(); i++)
    {
        elems_[i].id = i;
        elems_[i].neigh_ptrs.resize(4, nullptr);
        elems_[i].type = _4NodeTetrahedron;

        int_tag[i][2] = 3; // 剖分组别为三维

        elems_[i].vert_ptrs.resize(4);
        elems_[i].vert_ptrs[0] = nodes_.get(tet_3d[i].vert[0]->id);
        elems_[i].vert_ptrs[1] = nodes_.get(tet_3d[i].vert[1]->id);
        elems_[i].vert_ptrs[2] = nodes_.get(tet_3d[i].vert[2]->id);
        elems_[i].vert_ptrs[3] = nodes_.get(tet_3d[i].vert[3]->id);
    }

    int offset =  tet_3d.size();
    for (size_t i = 0; i < tri_3d.size(); i++)
    {;
        elems_[i + offset].id = i + offset;
        elems_[i + offset].neigh_ptrs.resize(3, nullptr);
        elems_[i + offset].type = _3NodeTriangle;

        int_tag[i + offset][2] = 2;

        elems_[i + offset].vert_ptrs.resize(3);
        elems_[i + offset].vert_ptrs[0] = nodes_.get(tri_3d[i].vert[0]->id);
        elems_[i + offset].vert_ptrs[1] = nodes_.get(tri_3d[i].vert[1]->id);
        elems_[i + offset].vert_ptrs[2] = nodes_.get(tri_3d[i].vert[2]->id);
    }

    if (!tet_tag.empty())
    {
        for (size_t i = 0; i < tet_3d.size(); i++)
        {
            int_tag[i][1] = tet_tag[i];
        }
    }

    if (!tri_tag.empty())
    {
        for (size_t i = 0; i < tri_3d.size(); i++)
        {
            int_tag[i + offset][1] = tri_tag[i];
        }
    }
    
    // 整理单元体组
    meshio_element_group tmp_group;
    tmp_group.type = elems_[0].type;
    tmp_group.phys_group = int_tag[0][0];
    tmp_group.geom_group = int_tag[0][1];
    tmp_group.part_group = int_tag[0][2];
    tmp_group.elem_ptrs.push_back(elems_.get(0));
    groups_.push_back(tmp_group);

    bool not_found;
    for (size_t i = 1; i < elems_.size(); i++)
    {
        not_found = true;
        for (size_t g = 0; g < groups_.size(); g++)
        {
            if (groups_[g].type == elems_[i].type &&
                groups_[g].phys_group == int_tag[i][0] &&
                groups_[g].geom_group == int_tag[i][1] &&
                groups_[g].part_group == int_tag[i][2])
            {
                groups_[g].elem_ptrs.push_back(elems_.get(i));
                not_found = false;
                break;
            }
        }
        
        if (not_found)
        {
            tmp_group.elem_ptrs.clear();
            tmp_group.type = elems_[i].type;
            tmp_group.phys_group = int_tag[i][0];
            tmp_group.geom_group = int_tag[i][1];
            tmp_group.part_group = int_tag[i][2];
            tmp_group.elem_ptrs.push_back(elems_.get(i));
            groups_.push_back(tmp_group);
        }
    }

    std::string n_file = filename + ".neigh";
    if (!access(n_file.c_str(), R_OK))
    {
        read_Tetgen_neighbor(filename, tet_3d, is_packed);

        for (size_t i = 0; i < tet_3d.size(); i++)
        {
            for (size_t j = 0; j < 4; j++)
            {
                if (tet_3d[i].neigh[j] != nullptr) elems_[i].neigh_ptrs[j] = elems_.get(tet_3d[i].neigh[j]->id);
            }
        }
    }

    // 最后将所有的组别和单元体设置为可用
    valid_group_size_ = groups_.size();
    for (size_t g = 0; g < groups_.size(); g++)
    {
        groups_[g].enabled = true;
        groups_[g].enable_elements();
    }
    initialized_ = true;
    return;
}

void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_packed)
{
    if (initialized_ == true)
    {
        throw std::runtime_error("[gctl::mesh_io::read_gmsh_v2_ascii()] The mesh is already initialized. Run clear() first.");
    }

    std::ifstream infile;
    open_infile(infile, filename, ".msh");

    std::string tmp_str, tmp_str2;
    std::stringstream tmp_ss;

    // 检查mesh文件版本
    while(getline(infile,tmp_str))
    {
        if (tmp_str == "$MeshFormat")
        {
            getline(infile, tmp_str);
            if (tmp_str != "2.2 0 8")
            {
                throw std::runtime_error("[gctl::mesh_io::read_gmsh_v2_ascii()] Unsupported mesh file format.");
            }
            break;
        }
    }

    //读入模型空间顶点集和物理分组
    size_t p_size = 0;
    size_t n_size = 0;
    array<gmsh_physical_group> phys;

    while(getline(infile,tmp_str))
    {
        if (tmp_str == "$PhysicalNames")
        {
            getline(infile, tmp_str);
            gctl::str2ss(tmp_str, tmp_ss);
            tmp_ss >> p_size;

            phys.resize(p_size);
            for (size_t i = 0; i < p_size; i++)
            {
                getline(infile, tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> phys[i].dim_tag >> phys[i].phys_tag >> tmp_str2;
                // 去掉两端的双引号
                replace_all(phys[i].name, tmp_str2, "\"", "");
            }
        }

        if (tmp_str == "$Nodes")
        {
            getline(infile, tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> n_size; //第一个数为顶点个数
            nodes_.resize(n_size);
            valid_node_size_ = n_size;

            for (int i = 0; i < n_size; i++)
            {
                getline(infile, tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> nodes_[i].id >> nodes_[i].x >> nodes_[i].y >> nodes_[i].z;
                if (is_packed == NotPacked) nodes_[i].id -= 1;
            }
            break;
        }
    }

    // 读入模型空间元素集
    int i_size, type_code, attri_num, vt_idx;
    array<array<int> > file_itag;
    while(getline(infile,tmp_str))
    {
        if (tmp_str == "$Elements")
        {
            getline(infile,tmp_str);
            gctl::str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size; //第一个数为元素个数

            valid_elem_size_ = i_size;
            elems_.resize(i_size);
            file_itag.resize(valid_elem_size_);

            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile,tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> elems_[i].id >> type_code >> attri_num;
                if (is_packed == NotPacked) elems_[i].id -= 1;

                elems_[i].type = elem_gmsh_type(type_code);
                elems_[i].vert_ptrs.resize(elem_size(elems_[i].type));

                // we get at least three tags. see below for tag types
                // default tags will be assgined to DEFAULT_INVALID_TAG
                file_itag[i].resize(GCTL_MAX(3, attri_num), DEFAULT_INVALID_TAG);
                for (size_t a = 0; a < attri_num; a++)
                {
                    tmp_ss >> file_itag[i][a];
                }

                for (size_t v = 0; v < elems_[i].vert_ptrs.size(); v++)
                {
                    tmp_ss >> vt_idx;
                    if (is_packed == Packed) elems_[i].vert_ptrs[v] = nodes_.get(vt_idx);
                    else elems_[i].vert_ptrs[v] = nodes_.get(vt_idx - 1);
                }
            }
            break;
        }
    }

    // 读入数据模块
    meshio_data tmp_data;
    while(getline(infile,tmp_str))
    {
        if (!infile.good()) break;

        if (tmp_str == "$NodeData" || tmp_str == "$ElementData")
        {
            tmp_data.clear();
            tmp_data.enabled = true;

            if (tmp_str == "$NodeData") tmp_data.d_type = NodeData;
            if (tmp_str == "$ElementData") tmp_data.d_type = ElemData;

            getline(infile,tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size;
            tmp_data.str_tag.resize(i_size);
            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile, tmp_str);
                replace_all(tmp_str2, tmp_str, "\"", "");
                tmp_data.str_tag[i] = tmp_str2;
            }

            getline(infile,tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size;
            tmp_data.real_tag.resize(i_size);
            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile,tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> tmp_data.real_tag[i];
            }

            getline(infile,tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size;
            tmp_data.int_tag.resize(i_size);
            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile,tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> tmp_data.int_tag[i];
            }

            if (tmp_data.d_type == NodeData)
            {
                tmp_data.vert_ptrs.resize(tmp_data.int_tag.back());
                tmp_data.val.resize(tmp_data.int_tag.back());
                for (size_t i = 0; i < tmp_data.val.size(); i++)
                {
                    getline(infile,tmp_str);
                    str2ss(tmp_str, tmp_ss);
                    tmp_ss >> vt_idx >> tmp_data.val[i];
                    if (is_packed == Packed) tmp_data.vert_ptrs[i] = nodes_.get(vt_idx);
                    else tmp_data.vert_ptrs[i] = nodes_.get(vt_idx - 1);
                }
            }
            else
            {
                tmp_data.elem_ptrs.resize(tmp_data.int_tag.back());
                tmp_data.val.resize(tmp_data.int_tag.back());
                for (size_t i = 0; i < tmp_data.val.size(); i++)
                {
                    getline(infile,tmp_str);
                    str2ss(tmp_str, tmp_ss);
                    tmp_ss >> vt_idx >> tmp_data.val[i];
                    if (is_packed == Packed) tmp_data.elem_ptrs[i] = elems_.get(vt_idx);
                    else tmp_data.elem_ptrs[i] = elems_.get(vt_idx - 1);
                }
            }

            datas_.push_back(tmp_data);
            continue;
        }
    }
    infile.close();

    // 整理单元体组
    meshio_element_group tmp_group;
    tmp_group.type = elems_[0].type;
    tmp_group.phys_group = file_itag[0][0];
    tmp_group.geom_group = file_itag[0][1];
    tmp_group.part_group = file_itag[0][2];
    tmp_group.elem_ptrs.push_back(elems_.get(0));
    groups_.push_back(tmp_group);

    // 元素类型与三个标记都一致的元素将被分到同一组
    bool not_found;
    for (size_t i = 1; i < elems_.size(); i++)
    {
        not_found = true;
        for (size_t g = 0; g < groups_.size(); g++)
        {
            if (groups_[g].type == elems_[i].type &&
                groups_[g].phys_group == file_itag[i][0] &&
                groups_[g].geom_group == file_itag[i][1] &&
                groups_[g].part_group == file_itag[i][2])
            {
                groups_[g].elem_ptrs.push_back(elems_.get(i));
                not_found = false;
                break;
            }
        }
        
        if (not_found)
        {
            tmp_group.elem_ptrs.clear();
            tmp_group.type = elems_[i].type;
            tmp_group.phys_group = file_itag[i][0]; // 物理组
            tmp_group.geom_group = file_itag[i][1]; // 几何组
            tmp_group.part_group = file_itag[i][2]; // 剖分组（一般以元素的维度区分）
            tmp_group.elem_ptrs.push_back(elems_.get(i));
            groups_.push_back(tmp_group);
        }
    }

    if (!phys.empty())
    {
        // 遍历所有元素组 按物理组为标准为元素组命名 并以将元素维度赋值为剖分组
        for (size_t g = 0; g < groups_.size(); g++)
        {
            for (size_t p = 0; p < phys.size(); p++)
            {
                if (phys[p].phys_tag == groups_[g].phys_group)
                {
                    groups_[g].name = phys[p].name;
                    groups_[g].part_group = phys[p].dim_tag;
                    break;
                }
            }
        }
    }

    // 最后将所有的组别和单元体设置为可用
    valid_group_size_ = groups_.size();
    for (size_t g = 0; g < groups_.size(); g++)
    {
        groups_[g].enabled = true;
        groups_[g].enable_elements();
    }
    initialized_ = true;
    return;
}

void gctl::mesh_io::save_gmsh_v2_ascii(std::string filename, index_packed_e is_packed)
{
    if (initialized_ == false)
    {
        throw std::runtime_error("The object is not initialized. From gctl::mesh_io::save_mesh(...)");
    }

    std::ofstream outfile;
    open_outfile(outfile, filename, ".msh");

    outfile << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n";

    std::vector<gmsh_physical_group> gmsh_groups;
    get_gmsh_physical_groups(gmsh_groups);

    outfile << "$PhysicalNames\n" << gmsh_groups.size() << "\n";
    for (size_t i = 0; i < gmsh_groups.size(); i++)
    {
        if (gmsh_groups[i].phys_tag != DEFAULT_INVALID_TAG) outfile << gmsh_groups[i].dim_tag << " " << gmsh_groups[i].phys_tag << " \"" << gmsh_groups[i].name << "\"\n";
        else outfile << gmsh_groups[i].dim_tag << " 0 \"" << gmsh_groups[i].name << "\"\n";
    }
    outfile << "$EndPhysicalNames\n";

    outfile << "$Nodes\n" << valid_node_size_ << std::endl;
    for (size_t i = 0; i < nodes_.size(); i++)
    {
        if (nodes_[i].id == DEFAULT_INVALID_TAG) continue;
    
        if (is_packed == NotPacked) outfile << nodes_[i].id + 1 << " ";
        else outfile << nodes_[i].id << " ";

        outfile << std::setprecision(16) << nodes_[i].x << " " << nodes_[i].y << " " << nodes_[i].z << std::endl;
    }
    outfile<<"$EndNodes\n";

    size_t a_size;
    outfile << "$Elements\n" << valid_elem_size_ << std::endl;
    for (size_t g = 0; g < groups_.size(); g++)
    {
        if (groups_[g].enabled)
        {
            for (size_t e = 0; e < groups_[g].elem_ptrs.size(); e++)
            {
                if (is_packed == NotPacked) outfile << groups_[g].elem_ptrs[e]->id + 1 << " " << elem_gmsh_code(groups_[g].type) << " ";
                else outfile << groups_[g].elem_ptrs[e]->id << " " << elem_gmsh_code(groups_[g].type) << " ";

                a_size = 1;
                if (groups_[g].geom_group != DEFAULT_INVALID_TAG) a_size++;
                if (groups_[g].part_group != DEFAULT_INVALID_TAG) a_size++;

                outfile << a_size;
                if (groups_[g].phys_group != DEFAULT_INVALID_TAG) outfile << " " << groups_[g].phys_group;
                else outfile << " 0";

                if (groups_[g].geom_group != DEFAULT_INVALID_TAG) outfile << " " << groups_[g].geom_group;
                if (groups_[g].part_group != DEFAULT_INVALID_TAG) outfile << " " << groups_[g].part_group;

                for (size_t v = 0; v < groups_[g].elem_ptrs[e]->vert_ptrs.size(); v++)
                {
                    if (is_packed == NotPacked) outfile << " " << groups_[g].elem_ptrs[e]->vert_ptrs[v]->id + 1;
                    else outfile << " " << groups_[g].elem_ptrs[e]->vert_ptrs[v]->id;
                }
                outfile << std::endl;
            }
        }
    }

    outfile << "$EndElements\n";

    if (!datas_.empty())
    {
        size_t d_size;
        for (size_t i = 0; i < datas_.size(); i++)
        {
            if (!datas_[i].enabled) continue;

            if (datas_[i].d_type == NodeData)
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].vert_ptrs[n]->id != DEFAULT_INVALID_TAG) d_size++;
                }
                
                if (d_size)
                {
                    outfile << "$NodeData\n";
                    outfile << datas_[i].str_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].str_tag.size(); a++)
                    {
                        outfile << "\"" << datas_[i].str_tag[a] << "\"\n";
                    }

                    outfile << datas_[i].real_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].real_tag.size(); a++)
                    {
                        outfile << datas_[i].real_tag[a] << "\n";
                    }

                    outfile << datas_[i].int_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].int_tag.size() - 1; a++)
                    {
                        outfile << datas_[i].int_tag[a] << "\n";
                    }
                    outfile << d_size << "\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].vert_ptrs[a]->id != DEFAULT_INVALID_TAG && is_packed == Packed) outfile << datas_[i].vert_ptrs[a]->id << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                        else if (datas_[i].vert_ptrs[a]->id != DEFAULT_INVALID_TAG) outfile << datas_[i].vert_ptrs[a]->id + 1 << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                    }
                    outfile << "$EndNodeData\n";
                }
            }
            else
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].elem_ptrs[n]->enabled) d_size++;
                }
                
                if (d_size)
                {
                    outfile << "$ElementData\n";
                    outfile << datas_[i].str_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].str_tag.size(); a++)
                    {
                        outfile << "\"" << datas_[i].str_tag[a] << "\"\n";
                    }

                    outfile << datas_[i].real_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].real_tag.size(); a++)
                    {
                        outfile << datas_[i].real_tag[a] << "\n";
                    }

                    outfile << datas_[i].int_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].int_tag.size() - 1; a++)
                    {
                        outfile << datas_[i].int_tag[a] << "\n";
                    }
                    outfile << d_size << "\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].elem_ptrs[a]->enabled && is_packed == Packed) outfile << datas_[i].elem_ptrs[a]->id << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                        else if (datas_[i].elem_ptrs[a]->enabled) outfile << datas_[i].elem_ptrs[a]->id + 1 << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                    }
                    outfile << "$EndElementData\n";
                }
            }
        }
    }

    outfile.close();
    return;
}

void gctl::mesh_io::save_vtk_legacy_ascii(std::string filename)
{
    if (initialized_ == false)
    {
        throw std::runtime_error("The object is not initialized. From gctl::mesh_io::export_vtk(...)");
    }

    std::ofstream outfile;
    open_outfile(outfile, filename, ".vtk");

    outfile << "# vtk DataFile Version 2.0\nGenerated by gctl::mesh_io\nASCII\nDATASET UNSTRUCTURED_GRID\nPOINTS " << valid_node_size_ <<" FLOAT\n";
	for (int i = 0; i < nodes_.size(); i++)
	{
        if (nodes_[i].id == DEFAULT_INVALID_TAG) continue;

		outfile << std::setprecision(16) << nodes_[i].x << " " << nodes_[i].y << " " << nodes_[i].z << std::endl;
	}

	//计算一下CELLS的总长 注意类型名算一个 所以要加1
    size_t totalCellNum = 0;
	for (int i = 0; i < elems_.size(); i++)
	{
        if (elems_[i].enabled) totalCellNum += elem_size(elems_[i].type) + 1;
	}
	outfile << "CELLS " << valid_elem_size_ << " " << totalCellNum << std::endl;

	for (int i = 0; i < elems_.size(); i++)
	{
        if (!elems_[i].enabled) continue;

		outfile << elem_size(elems_[i].type) << " ";
		for (int j = 0; j < elems_[i].vert_ptrs.size(); j++)
		{
			outfile << elems_[i].vert_ptrs[j]->id << " ";
		}
		outfile << std::endl;
	}

	outfile << "CELL_TYPES " << valid_elem_size_ << std::endl;
	for (int i = 0; i < elems_.size(); i++)
	{
		if (elems_[i].enabled) outfile << elem_vtk_code(elems_[i].type) << "\n";
	}

    bool nd_head = true, ed_head = true;
    size_t d_size;
    std::string temp_dataName;
    if (!datas_.empty())
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            if (!datas_[i].enabled) continue;

            if (datas_[i].d_type == NodeData)
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].vert_ptrs[n]->id != DEFAULT_INVALID_TAG) d_size++;
                }
                
                if (d_size && d_size == valid_node_size_)
                {
                    if (nd_head) {outfile<<"POINT_DATA "<< d_size << std::endl; nd_head = false; ed_head = true;}

                    gctl::replace_all(temp_dataName, datas_[i].str_tag.front()," ","_");
                    outfile << "SCALARS " << temp_dataName << " FLOAT\nLOOKUP_TABLE default\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].vert_ptrs[a]->id != DEFAULT_INVALID_TAG) outfile << std::setprecision(16) <<  datas_[i].val[a] << "\n";
                    }
                }
            }
            else
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].elem_ptrs[n]->enabled) d_size++;
                }
                
                if (d_size && d_size == valid_elem_size_)
                {
                    if (ed_head) {outfile<<"CELL_DATA "<< d_size << std::endl; ed_head = false; nd_head = true;}

                    gctl::replace_all(temp_dataName, datas_[i].str_tag.front()," ","_");
                    outfile << "SCALARS " << temp_dataName << " FLOAT\nLOOKUP_TABLE default\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].elem_ptrs[a]->enabled) outfile << std::setprecision(16) <<  datas_[i].val[a] << "\n";
                    }
                }
            }
        }
    }

    outfile.close();
    return;
}

void gctl::mesh_io::save_data_to_xyz(std::string filename, std::string dataname, coordinate_system_e out_coor, double refr, double refR)
{
    int tmp_size;
    std::string tmp_name;
    point3ds ps;
    point3dc pc;

    std::ofstream ofile;
    for (size_t i = 0; i < datas_.size(); i++)
    {
        if (dataname != "null" && (datas_[i].str_tag.front() != dataname || !datas_[i].enabled)) continue;

        gctl::replace_all(tmp_name, datas_[i].str_tag.front()," ","_");
        gctl::replace_all(tmp_name, tmp_name,"/","");
        open_outfile(ofile, filename + "_" + tmp_name, ".txt");

        if (datas_[i].d_type == NodeData)
        {
            for (size_t n = 0; n < datas_[i].vert_ptrs.size(); n++)
            {
                if (datas_[i].vert_ptrs[n]->id != DEFAULT_INVALID_TAG)
                {
                    if (out_coor == Spherical)
                    {
                        ps = datas_[i].vert_ptrs[n]->c2s();
                        ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*M_PI/180.0);

                        ofile << ps.lon << " " << ps.lat << " " << ps.rad << " " << datas_[i].val[n] << "\n";
                    }
                    else  ofile << pc.x << " " << pc.y << " " << pc.z << " " << datas_[i].val[n] << "\n";
                }
            }
            
        }
        else
        {
            for (size_t e = 0; e < datas_[i].elem_ptrs.size(); e++)
            {
                if (datas_[i].elem_ptrs[e]->enabled)
                {
                    pc = point3dc(0.0, 0.0, 0.0);
                    tmp_size = datas_[i].elem_ptrs[e]->vert_ptrs.size();
                    for (size_t v = 0; v < datas_[i].elem_ptrs[e]->vert_ptrs.size(); v++)
                    {
                        pc.x += datas_[i].elem_ptrs[e]->vert_ptrs[v]->x;
                        pc.y += datas_[i].elem_ptrs[e]->vert_ptrs[v]->y;
                        pc.z += datas_[i].elem_ptrs[e]->vert_ptrs[v]->z;
                    }
                    pc.x /= tmp_size;
                    pc.y /= tmp_size;
                    pc.z /= tmp_size;

                    if (out_coor == Spherical)
                    {
                        ps = pc.c2s();
                        ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*M_PI/180.0);

                        ofile << ps.lon << " " << ps.lat << " " << ps.rad << " " << datas_[i].val[e] << "\n";
                    }
                    else ofile << pc.x << " " << pc.y << " " << pc.z << " " << datas_[i].val[e] << "\n";
                }
            }
        }
        ofile.close();
    }
    return;
}

std::string gctl::mesh_io::elem_name(element_type_enum e_type)
{
    return elem_type2name_[e_type];
}

int gctl::mesh_io::elem_gmsh_code(element_type_enum e_type)
{
    return elem_type2gmshcode_[e_type];
}

int gctl::mesh_io::elem_vtk_code(element_type_enum e_type)
{
    return elem_type2vtkcode_[e_type];
}

int gctl::mesh_io::elem_size(element_type_enum e_type)
{
    return elem_type2size_[e_type];
}

gctl::element_type_enum gctl::mesh_io::elem_gmsh_type(int code)
{
    return elem_gmshcode2type_[code];
}

gctl::element_type_enum gctl::mesh_io::elem_vtk_type(int code)
{
    return elem_vtkcode2type_[code];
}

void gctl::mesh_io::update_indexing()
{
    // 重置所有顶点的索引
    for (size_t n = 0; n < nodes_.size(); n++)
    {
        nodes_[n].id = DEFAULT_INVALID_TAG;
    }
    
    // 确定有效单元个数并初始化对应的顶点索引为1
    int new_id = 0;
    for (size_t e = 0; e < elems_.size(); e++)
    {
        if (elems_[e].enabled)
        {
            elems_[e].id = new_id; new_id++;
            for (size_t n = 0; n < elems_[e].vert_ptrs.size(); n++)
            {
                elems_[e].vert_ptrs[n]->id = 1;
            }
        }
    }
    valid_elem_size_ = new_id;

    // 对所有有效顶点赋新的索引 同时确定有效顶点的个数
    new_id = 0;
    for (size_t n = 0; n < nodes_.size(); n++)
    {
        if (nodes_[n].id == 1)
        {
            nodes_[n].id = new_id; new_id++;
        }
    }
    valid_node_size_ = new_id;

    // 确定有效单元组的个数
    valid_group_size_ = 0;
    for (size_t p = 0; p < groups_.size(); p++)
    {
        if (groups_[p].enabled) valid_group_size_++;
    }
    return;
}

void gctl::mesh_io::sort_groups()
{
    // 拷贝到临时组
    std::vector<meshio_element_group> tmp_groups = groups_;

    // 清空对象
    for (size_t i = 0; i < groups_.size(); i++)
    {
        groups_[i].elem_ptrs.clear();
    }
    groups_.clear();
    
    // 整理单元体组
    bool not_found;
    meshio_element_group tmp_group;
    for (size_t i = 0; i < tmp_groups.size(); i++)
    {
        tmp_group = tmp_groups[i];

        not_found = true;
        for (size_t g = 0; g < groups_.size(); g++)
        {
            // 注意名称不是分组的标准 多个组可以共用一个名称
            if (groups_[g].type == tmp_group.type &&
                groups_[g].phys_group == tmp_group.phys_group &&
                groups_[g].geom_group == tmp_group.geom_group &&
                groups_[g].part_group == tmp_group.part_group)
            {
                for (size_t e = 0; e < tmp_group.elem_ptrs.size(); e++)
                {
                    groups_[g].elem_ptrs.push_back(tmp_group.elem_ptrs[e]);
                }

                not_found = false;
            }
            
        }

        if (not_found)
        {
            groups_.push_back(tmp_group);
        }

        tmp_group.elem_ptrs.clear();
    }

    valid_group_size_ = 0;
    for (size_t i = 0; i < groups_.size(); i++)
    {
        if (groups_[i].enabled) valid_group_size_++;
    }
    return;
}