#include "tetgenext.h"
#if defined(_WIN32)
#include <Windows.h>
#endif


namespace TETGEN_EXT {
    void DebugBDF(std::ostream& out, size_t nnode, size_t nelem,
        std::function<bool(size_t i, int32_t &nid, double *pos)> nodeFun,
        std::function<bool(size_t i, int32_t &elemId, int32_t &group, int32_t &, int32_t *nodes)> elemFun)
    {
        auto  double2Char16 = [&](double val, char* buffer)
        {
            sprintf(buffer, "%16.16lf", val);
            buffer[15] = '\0';
            return buffer;
        };

        auto  BdfPrintNodeLong = [&](int id, double x, double y, double z, char* strBuffer)
        {
            char bufx[32], bufy[32], bufz[32];
            sprintf(strBuffer, "%-16s%8d%-16s%-16s%-16s%-1s\n%-8s%-16s\n", "GRID*", id, "",
                double2Char16(x, bufx), double2Char16(y, bufy), "*", "*", double2Char16(z, bufz));
        };


        auto WriteElementBdf = [&](std::ostream& out, int id, int tp, int gp, const int* nodes0, int nid_offset)
        {
            out.setf(std::ios::left);
            int num = tp % 100;
            int nodes[20];
            if (gp == 0)
                gp = tp;
            memcpy(nodes, nodes0, num * sizeof(int));
            if (tp == 203)
            {
                out << "CTRIA3  " << std::setw(8) << id << std::setw(8) << gp;
            }
            else if (tp == 204)
            {
                out << "CQUAD4  " << std::setw(8) << id << std::setw(8) << gp;
            }
            else if (tp == 304)
            {
                out << "CTETRA  " << std::setw(8) << id << std::setw(8) << gp;
            }
            else if (tp == 308)
            {
                out << "CHEXA   " << std::setw(8) << id << std::setw(8) << gp;
            }
            else if (tp == 306)
            {
                out << "CPENTA  " << std::setw(8) << id << std::setw(8) << gp;
            }
            else if (tp == 305)
            {
                //out << "CPYRA  " << std::setw(8) << id << std::setw(8) << gp;
                out << "CHEXA   " << std::setw(8) << id << std::setw(8) << gp;
                tp = 308;
                nodes[7] = nodes[6] = nodes[5] = nodes[4];
                num = 8;
            }
            else if (tp == 102)
            {
                if (gp == 0)
                    gp = 1;
                // out<<  "PLOTEL  "<<std::setw(8)<<id;
                out << "CBAR    " << std::setw(8) << id << std::setw(8) << gp;
                // out<<  "PLOTEL  "<<","<<id;
            }



            if (num > 0)
            {
                int ncount = 0;
                int jfield = 0;
                // first line
                for (jfield = 4; jfield <= 9; jfield++)
                {
                    ncount++;
                    if (ncount > num) break;
                    out << std::setw(8) << nodes[ncount - 1] + nid_offset;
                }
                out << std::endl;

                // following lines        
                while (ncount < num)
                {
                    out << std::setw(8) << "+";
                    for (jfield = 2; jfield <= 9; jfield++)
                    {
                        ncount++;
                        if (ncount > num) break;
                        out << std::setw(8) << nodes[ncount - 1] + nid_offset;
                    }
                    out << std::endl;
                }

            } // if (num>0)     

        };


        auto  BdfPrintNodeLong2 = [&](int id, double x, double y, double z, std::ostream& out)
        {
            /*out << "GRID," << id << ",,"
            << boost::str(boost::format("%1$g") % x) << ","
            << boost::str(boost::format("%1$g") % y) << ","
            << boost::str(boost::format("%1$g") % z) << std::endl;*/
            char strBuffer[1024];
            BdfPrintNodeLong(id, x, y, z, strBuffer);
            out << strBuffer;
        };

        int32_t nid, elemId, group, type, nodes[20];
        double pos[3];
        for (size_t i = 0; i < nnode; i++) {
            if (!nodeFun(i, nid, pos)) {
                continue;
            }
            BdfPrintNodeLong2(nid, pos[0], pos[1], pos[2], out);
        }

        for (size_t i = 0; i < nelem; i++) {
            if (!elemFun(i, elemId, group, type, nodes)) {
                continue;
            }
            WriteElementBdf(out, elemId, type, group, nodes, 0);
        }
    }

    void DebugBDF(const std::string & fn, size_t nnode, size_t nelem,
        std::function<bool(size_t i, int32_t &, double *)> nodeFun,
        std::function<bool(size_t i, int32_t &, int32_t &, int32_t &, int32_t *)> elemFun)
    {
        std::ofstream out(fn);
        DebugBDF(out, nnode, nelem, nodeFun, elemFun);
    }


    void DebugBdf(const tetgenio& out, const char* fn, const std::vector<int>* flagTet, bool asGroupId = false)
    {
        DebugBDF(fn, out.numberofpoints, out.numberoftetrahedra + out.numberoffacets, [&out](size_t i, int32_t& nid, double* pos) {
            nid = (int32_t)i + 1;
            memcpy(pos, out.pointlist + i * 3, 3 * sizeof(double));
            return true;
            }, [&](size_t i, int32_t& elemId, int32_t& group, int32_t&type, int32_t* nodes) {
                elemId = (int32_t)i + 1;
                if (i < out.numberoftetrahedra)
                {
                    type = 304;
                    group = type;
                    auto tet = out.tetrahedronlist + i * 4;
                    nodes[0] = tet[0] + 1;
                    nodes[1] = tet[1] + 1;
                    nodes[2] = tet[2] + 1;
                    nodes[3] = tet[3] + 1;
                    if (flagTet)
                    {
                        if (asGroupId)
                        {
                            group = flagTet->at(i);
                        }
                        else if (flagTet->at(i) < 0) {
                            group *= 10;
                        }
                    }
                }
                else
                {
                    type = 203;
                    group = type;
                    i -= out.numberoftetrahedra;
                    const auto& p = out.facetlist[i].polygonlist[0];
                    nodes[0] = p.vertexlist[0] + 1;
                    nodes[1] = p.vertexlist[1] + 1;
                    nodes[2] = p.vertexlist[2] + 1;
                }
                return true;
            });
    }

    static double GetMaxSize(Json::Value& meshCtrl)
    {
        if (meshCtrl.isMember("meshSizeMax")) {
            double dfmax;
            std::stringstream ss2;
            ss2 << meshCtrl["meshSizeMax"].asString();
            ss2 >> dfmax;
            return dfmax;
        }
        return -1.0;
    }

    static std::string Convert2TetGenCommand(Json::Value& meshCtrl, unsigned nInputTet = 0)
    {
        /*
        -p	Tetrahedralizes a piecewise linear complex (PLC).
        -Y	Preserves the input surface mesh (does not modify it).
        -r	Reconstructs a previously generated mesh.
        -q	Refines mesh (to improve mesh quality).
        -R	Mesh coarsening (to reduce the mesh elements).
        -A	Assigns attributes to tetrahedra in different regions.
        -a	Applies a maximum tetrahedron volume constraint.
        -m	Applies a mesh sizing function.
        -i	Inserts a list of additional points.
        -O	Specifies the level of mesh optimization.
        -S	Specifies maximum number of added points.
        -T	Sets a tolerance for coplanar test (default 10−8).
        -X	Suppresses use of exact arithmetic.
        -M	No merge of coplanar facets or very close vertices.
        -w	Generates weighted Delaunay (regular) triangulation.
        -c	Retains the convex hull of the PLC.
        -d	Detects self-intersections of facets of the PLC.
        -z	Numbers all output items starting from zero.
        -f	Outputs all faces to .face file.
        -e	Outputs all edges to .edge file.
        -n	Outputs tetrahedra neighbors to .neigh file.
        -v	Outputs Voronoi diagram to files.
        -g	Outputs mesh to .mesh file for viewing by Medit.
        -k	Outputs mesh to .vtk file for viewing by Paraview.
        -J	No jettison of unused vertices from output .node file.
        -B	Suppresses output of boundary information.
        -N	Suppresses output of .node file.
        -E	Suppresses output of .ele file.
        -F	Suppresses output of .face and .edge file.
        -I	Suppresses mesh iteration numbers.
        -C	Checks the consistency of the final mesh.
        -Q	Quiet: No terminal output except errors.
        -V	Verbose: Detailed information, more terminal output.
        -h	Help: A brief instruction for using TetGen.
        */
        std::stringstream ss;
        if (meshCtrl.isMember("initialDelaunaySimple")) {
            ss << "pn";
            return ss.str();
        }
        ss << "Ypni";
        //ss << "pni"; //for tetrefine.
        if (meshCtrl.empty()) {
            ss << "q1.2CC";
            return ss.str();
        }

        double dfmax = GetMaxSize(meshCtrl);
        if (dfmax > 1.e-10) {
            dfmax *= 1.4;
            double volMax = dfmax * dfmax*dfmax*0.118;
            ss << "a" << volMax;
        }

        if (meshCtrl.isMember("convexHull")) {
            ss << "c";
            return ss.str();
        }
        else if (meshCtrl.isMember("initialDelaunay")) {
            Json::Value optCtrl = meshCtrl["initialDelaunay"];
            if (optCtrl.isMember("removeSteiners")) {
                ss << "R";
            }
            else if (optCtrl.isMember("removeSteiners4BL"))
            {
                ss << "R";
                ss << "Y///1";
            }
        }
        else if (meshCtrl.isMember("optimize")) {
            if (nInputTet > 0)
                ss << "r";
            Json::Value optCtrl = meshCtrl["optimize"];
            if (optCtrl.isMember("operationFlag")) {
                int flag = optCtrl["operationFlag"].asInt();
                if (flag > 0) {
                    ss << "O2/" << flag;
                }
            }
            if (optCtrl.isMember("dihedralAngle")) {
                double val = optCtrl["dihedralAngle"].asDouble();
                if (val > 90) {
                    ss << "o/" << val;
                }
            }

            double reR = 1.414, daM = 0;
            if (optCtrl.isMember("RERatio")) {
                double val = optCtrl["RERatio"].asDouble();
                if (val > 0) {
                    reR = val;
                }
            }
            if (optCtrl.isMember("dihedralAngleMin")) {
                double val = optCtrl["dihedralAngleMin"].asDouble();
                if (val < 90) {
                    daM = val;
                }
            }

            ss << "q" << reR << "/" << daM;
            return ss.str();
        }
        else {
            ss << "q1.2CC";
            return ss.str();
        }
        return ss.str();
    }

    bool ReadBdfAndWriteJson(const std::string &fn, const std::string & fnout)
    {
        using namespace std;
        auto getblock = [](const std::string& str, int which, int length, std::stringstream& sstr)
        {//创建输入流sstr
            std::string block;
            block = str.substr(8 * which, 8 * length);//substr()返回本字符串的一个子串，从index开始，长num个字符。
            sstr.str(block);
        };

        typedef struct Node_
        {
            size_t     id;
            double        xyz[3];
        }Node;

        typedef struct Elem_
        {
            size_t        id;
            int           type;
            int           nodes[8];
            int           geoId;
        }Cell;

        size_t Nnodes = 0, Ntris = 0, Nquads = 0, Ntets;
        vector<Node> nodes_in;
        vector<Cell> elems_in;
        vector<double> nodes;
        vector<int> tris;
        vector<int> quads;
        vector<int> tets;
        ifstream inf(fn.c_str());

        if (!inf)
        {
            std::cout << "input error in bdfrw" << std::endl;
            return false;
        }
        string  block;
        string buf;
        stringstream sstr;
        std::set<int> ids;
        while (inf.eof() == false)
        {
            sstr.str("");
            sstr.clear();
            getline(inf, buf);
            getblock(buf, 0, 1, sstr);
            block.clear();
            sstr >> block;
            /*    cout<<block<<endl;*/
            if (block == "GRID*")
            {

                Node newN;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 2, sstr);
                sstr >> newN.id;
                {
                    sstr.str("");
                    sstr.clear();
                    getblock(buf, 5, 2, sstr);
                    sstr >> newN.xyz[0];
                    if (sstr.eof() == false)
                    {
                        int a = 0;
                        sstr >> a;
                        if (a < 0)
                            newN.xyz[0] = std::pow(10, a)*newN.xyz[0];
                    }
                    sstr.str("");
                    sstr.clear();
                    getblock(buf, 7, 2, sstr);
                    sstr >> newN.xyz[1];
                    if (sstr.eof() == false)
                    {
                        int a = 0;
                        sstr >> a;
                        if (a < 0)
                            newN.xyz[1] = std::pow(10, a)*newN.xyz[1];
                    }
                    sstr.str("");
                    sstr.clear();
                    getline(inf, buf, '\n');
                    getblock(buf, 1, 2, sstr);
                    sstr >> newN.xyz[2];
                    if (sstr.eof() == false)
                    {
                        int a = 0;
                        sstr >> a;
                        if (a < 0)
                            newN.xyz[2] = std::pow(10, a)*newN.xyz[2];
                    }


                    nodes_in.push_back(newN);
                    Nnodes++;
                }

            }

            if (block == "CTRIA3*")
            {

                Cell newC;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 2, sstr);
                sstr >> newC.id;
                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 2, sstr);
                sstr >> newC.geoId;

                sstr.str("");
                sstr.clear();
                getblock(buf, 5, 2, sstr);
                sstr >> newC.nodes[0];
                sstr.str("");
                sstr.clear();
                getblock(buf, 7, 2, sstr);
                sstr >> newC.nodes[1];
                sstr.str("");
                sstr.clear();
                getline(inf, buf, '\n');
                getblock(buf, 1, 2, sstr);
                sstr >> newC.nodes[2];

                newC.type = 203;

                elems_in.push_back(newC);
                Ntris++;

            }

            if (block == "CQUAD4*")
            {

                Cell newC;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 2, sstr);
                sstr >> newC.id;
                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 2, sstr);
                sstr >> newC.geoId;

                sstr.str("");
                sstr.clear();
                getblock(buf, 5, 2, sstr);
                sstr >> newC.nodes[0];
                sstr.str("");
                sstr.clear();
                getblock(buf, 7, 2, sstr);
                sstr >> newC.nodes[1];
                sstr.str("");
                sstr.clear();
                getline(inf, buf, '\n');
                getblock(buf, 1, 2, sstr);
                sstr >> newC.nodes[2];
                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 2, sstr);
                sstr >> newC.nodes[3];

                newC.type = 204;

                elems_in.push_back(newC);
                Nquads++;

            }

            if (block == "CTETRA*")
            {
                Cell newC;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 2, sstr);
                sstr >> newC.id;
                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 2, sstr);
                sstr >> newC.geoId;

                sstr.str("");
                sstr.clear();
                getblock(buf, 5, 2, sstr);
                sstr >> newC.nodes[0];
                sstr.str("");
                sstr.clear();
                getblock(buf, 7, 2, sstr);
                sstr >> newC.nodes[1];
                sstr.str("");
                sstr.clear();
                getline(inf, buf, '\n');
                getblock(buf, 1, 2, sstr);
                sstr >> newC.nodes[2];
                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 2, sstr);
                sstr >> newC.nodes[3];

                newC.type = 304;

                elems_in.push_back(newC);
                Ntets++;

            }

            if (block == "GRID")
            {

                Node newN;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 1, sstr);
                sstr >> newN.id;
                {
                    sstr.str("");
                    sstr.clear();
                    getblock(buf, 3, 1, sstr);
                    sstr >> newN.xyz[0];
                    if (sstr.eof() == false)
                    {
                        int a = 0;
                        sstr >> a;
                        if (a < 0)
                            newN.xyz[0] = std::pow(10, a)*newN.xyz[0];
                    }
                    sstr.str("");
                    sstr.clear();
                    getblock(buf, 4, 1, sstr);
                    sstr >> newN.xyz[1];
                    if (sstr.eof() == false)
                    {
                        int a = 0;
                        sstr >> a;
                        if (a < 0)
                            newN.xyz[1] = std::pow(10, a)*newN.xyz[1];
                    }
                    sstr.str("");
                    sstr.clear();
                    getblock(buf, 5, 1, sstr);
                    sstr >> newN.xyz[2];
                    if (sstr.eof() == false)
                    {
                        int a = 0;
                        sstr >> a;
                        if (a < 0)
                            newN.xyz[2] = std::pow(10, a)*newN.xyz[2];
                    }

                    nodes_in.push_back(newN);
                    Nnodes++;
                }
            }

            if (block == "CTRIA3")
            {

                Cell newC;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 1, sstr);
                sstr >> newC.id;
                sstr.str("");
                sstr.clear();
                getblock(buf, 2, 1, sstr);
                sstr >> newC.geoId;

                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 1, sstr);
                sstr >> newC.nodes[0];
                sstr.str("");
                sstr.clear();
                getblock(buf, 4, 1, sstr);
                sstr >> newC.nodes[1];
                sstr.str("");
                sstr.clear();
                getblock(buf, 5, 1, sstr);
                sstr >> newC.nodes[2];

                newC.type = 203;

                elems_in.push_back(newC);
                Ntris++;


            }

            if (block == "CQUAD4")
            {

                Cell newC;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 1, sstr);
                sstr >> newC.id;
                sstr.str("");
                sstr.clear();
                getblock(buf, 2, 1, sstr);
                sstr >> newC.geoId;

                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 1, sstr);
                sstr >> newC.nodes[0];
                sstr.str("");
                sstr.clear();
                getblock(buf, 4, 1, sstr);
                sstr >> newC.nodes[1];
                sstr.str("");
                sstr.clear();
                getblock(buf, 5, 1, sstr);
                sstr >> newC.nodes[2];
                sstr.str("");
                sstr.clear();
                getblock(buf, 6, 1, sstr);
                sstr >> newC.nodes[3];

                newC.type = 204;

                elems_in.push_back(newC);
                Nquads++;

            }

            if (block == "CTETRA")
            {

                Cell newC;
                sstr.str("");
                sstr.clear();
                getblock(buf, 1, 1, sstr);
                sstr >> newC.id;
                sstr.str("");
                sstr.clear();
                getblock(buf, 2, 1, sstr);
                sstr >> newC.geoId;

                sstr.str("");
                sstr.clear();
                getblock(buf, 3, 1, sstr);
                sstr >> newC.nodes[0];
                sstr.str("");
                sstr.clear();
                getblock(buf, 4, 1, sstr);
                sstr >> newC.nodes[1];
                sstr.str("");
                sstr.clear();
                getblock(buf, 5, 1, sstr);
                sstr >> newC.nodes[2];
                sstr.str("");
                sstr.clear();
                getblock(buf, 6, 1, sstr);
                sstr >> newC.nodes[3];

                newC.type = 304;

                elems_in.push_back(newC);
                Ntets++;

            }

        }
        inf.close();


        int maxid = 1;
        for (size_t i = 0; i < nodes_in.size(); i++)
        {
            if (maxid < nodes_in[i].id)
            {
                maxid = (int)nodes_in[i].id;
            }
        }
        vector<int> idmap(maxid + 1, -1);
        int Nnode2 = 0;
        for (size_t i = 0; i < nodes_in.size(); ++i) {

            if (idmap[nodes_in[i].id] == -1)
            {
                for (int j = 0; j < 3; j++)
                    nodes.push_back(nodes_in[i].xyz[j]);
                idmap[nodes_in[i].id] = Nnode2;
                Nnode2++;
            }
        }

        int maxcellid = 1;
        for (size_t i = 0; i < elems_in.size(); i++)
        {
            if (maxcellid < elems_in[i].id)
            {
                maxcellid = elems_in[i].id;
            }
        }
        int Nelem2 = 0;
        vector<int> idcellmapMB(maxcellid + 1, -1);
        for (size_t i = 0; i < elems_in.size(); i++)
        {
            if (idcellmapMB[elems_in[i].id] == -1)
            {
                int geoId = elems_in[i].geoId;
                if (elems_in[i].type == 203)
                {
                    for (int j = 0; j < 3; j++)
                        tris.push_back(idmap[elems_in[i].nodes[j]]);
                    tris.push_back(geoId);
                }
                else if (elems_in[i].type == 204)
                {
                    for (int j = 0; j < 4; j++)
                        quads.push_back(idmap[elems_in[i].nodes[j]]);
                    quads.push_back(geoId);
                }
                else if (elems_in[i].type == 304)
                {
                    for (int j = 0; j < 4; j++)
                        tets.push_back(idmap[elems_in[i].nodes[j]]);
                    tets.push_back(geoId);
                }
                idcellmapMB[elems_in[i].id] = Nelem2;
                Nelem2++;
            }
        }

        //writeJson;
        if (1)
        {
            Json::Value jsonOut;
            if (!WriteData(jsonOut, "meshInput.points", nodes))
            {
                return false;
            }

            WriteData(jsonOut, "meshInput.triangles", tris);

            WriteData(jsonOut, "meshInput.quadrilaterals", quads);

            WriteData(jsonOut, "meshInput.tetrahedrons", tets);


            if (!SetToFile(jsonOut, fnout))
                return false;
        }

        cout << "Input BDF file successfully!" << endl;
        cout << "***************************************" << endl;
        cout << "Nnodes:         " << Nnode2 << endl;
        cout << "Ntris:          " << tris.size() / 4 << endl;
        cout << "Nquads:         " << quads.size() / 5 << endl;
        cout << "Ntets:         " << tets.size() / 5 << endl;
        cout << "***************************************" << endl;

        return true;

    }

#define PARSE_VA_MSG0(fmtO,lenO,fmtI) {va_list argptr; va_start(argptr,fmtI); vsprintf(fmtO,fmtI,argptr);va_end(argptr);}
#define PARSE_VA_MSG( fmtO,lenO,fmtI) char fmtO[lenO+1];PARSE_VA_MSG0(fmtO,lenO,fmtI);

    void DebugMessageBox(const char * mg, ...)
    {
#if defined(_WIN32) && defined(_DEBUG)
        if (!getenv("ENABLE_PAUSE_TETGEN"))
            return;
        char buffer[8192];
        PARSE_VA_MSG(buf, 8192, mg);
        DWORD pidwin = GetCurrentProcessId();
        sprintf(buffer, "%s (pid=%d)", buf, pidwin);
        MessageBoxA(GetDesktopWindow(), buffer, "Debuger", MB_OK);
#endif
    }
    struct RefineData
    {
        std::vector<float> mValue;
        double sizeData[4];
        int32_t  mSize[3];
        double mStart[3], mStep[3], mStepInv[3];
        double mDfMax;
        RefineData()
        {
            mSize[0] = mSize[1] = mSize[2] = 0.0;
        }
        size_t sizeYZ()const { return mSize[1] * mSize[2]; }
        size_t index(size_t x, size_t y, size_t z)const { return x * sizeYZ() + y * mSize[2] + z; }
        void getCellIndex(const REAL *pos, int32_t*idx) const
        {
            for (uint32_t i = 0; i < 3; i++) {

                idx[i] = (pos[i] - mStart[i]) * mStepInv[i];
                if (idx[i] >= (mSize[i] - 1))
                {
                    idx[i] = mSize[i] - 1;
                }
                if (idx[i] < 0) {
                    idx[i] = 0;
                }
            }
        }

        float valueI(int x, int y, int z)const { return mValue[index(x, y, z)]; }

        double operator()(const REAL* p)const
        {
            REAL t[3];
            int idx[3];
            for (int i = 0; i < 3; i++) {
                t[i] = (p[i] - mStart[i]) / mStep[0];
                if (t[i] < 2 || t[i] > (mSize[i] - 3)) {
                    this->getCellIndex(p, idx);
                    return mValue[index(idx[0], idx[1], idx[2])];
                }
            }

            int i = (int)(t[0]), j = (int)(t[1]), k = (int)(t[2]);
            double a = t[0] - i, b = t[1] - j, c = t[2] - k;
            double lVal = (1 - a) * (1 - b) * (1 - c) * valueI(i, j, k) + a * (1 - b) * (1 - c) * valueI(i + 1, j, k) +
                (1 - a) * (b) * (1 - c) * valueI(i, j + 1, k) + a * (b) * (1 - c) * valueI(i + 1, j + 1, k) +
                (1 - a) * (1 - b) * (c)*valueI(i, j, k + 1) + a * (1 - b) * (c)*valueI(i + 1, j, k + 1) +
                (1 - a) * (b) * (c)*valueI(i, j + 1, k + 1) + a * (b) * (c)*valueI(i + 1, j + 1, k + 1);
            if (lVal < sizeData[2])
                lVal = sizeData[2];
            if (lVal > sizeData[3])
                lVal = sizeData[3];

            return lVal;
        }
    };

    RefineData gRefineData;

    REAL Distance3D(const REAL* a, const REAL* b)
    {
        REAL dx = a[0] - b[0];
        REAL dy = a[1] - b[1];
        REAL dz = a[2] - b[2];
        return sqrt(dx * dx + dy * dy + dz * dz);
    }

    double  LenNormal(double len, double h0, double h1)
    {
        assert(h0 > 0 && h1 > 0);
        double diff = h0 - h1;
        double sum = h0 + h1;
        double det = sum * 0.5e-2;
        if (fabs(diff) < det)
        {
            return 2 * len / sum;
        }
        return len * log(h0 / h1) / diff;
    }

    bool IsTetraNeedRefine(REAL* v1, REAL* v2, REAL* v3, REAL* v4, REAL* el, REAL area)
    {
        static const  short TET_EDGE_VERTEX[6][2] = { {0,1},{1,2},{2,0},{0,3},{1,3},{2,3} };
        const REAL* pts[] = { v1 ,v2 ,v3 ,v4 };
        auto getMaxLength = [&](double& max)
        {
            max = 0;
            int whiMax = -1;
            for (int i = 0; i < 6; i++)
            {
                double dis = Distance3D(pts[TET_EDGE_VERTEX[i][0]], pts[TET_EDGE_VERTEX[i][1]]);
                if (dis > max)
                {
                    max = dis;
                    whiMax = i;
                }
            }
            return whiMax;
        };
        double  lMax;
        int whiMax = getMaxLength(lMax);

        if (gRefineData.mDfMax > 0 && lMax > gRefineData.mDfMax)
        {
            return true;
        }

        if (gRefineData.mValue.empty())
            return false;

        double dfs[2];
        dfs[0] = gRefineData(pts[TET_EDGE_VERTEX[whiMax][0]]);
        if (dfs[0] <= 0)
            return false;
        dfs[1] = gRefineData(pts[TET_EDGE_VERTEX[whiMax][1]]);
        if (dfs[1] <= 0)
            return false;

        double lenMax = LenNormal(lMax, dfs[0], dfs[1]);
        return lenMax > 1.5;
    }


    class  ScopedRedirect
    {
    public:
        ScopedRedirect(const char* fn = "logfile.log")
        {
            _savedStream = freopen(fn, "w", stdout);
        }
        ~ScopedRedirect()
        {
            fflush(_savedStream);
            _savedStream = freopen("CON", "w", stdout);
        }
    private:
        FILE *_savedStream;
    };


    

    int MeshGen2(const std::string& finput, const std::string& fout)
    {
        typedef ArrayHashProxy<DInt3Array, HashInt3> Int3HArray;

        double mDFMax = -1;


        try
        {
            auto  MarkTetra = [](const tetgenio& in, const tetgenio& out, DIntArray& tetFlag) {
                if (out.numberoftetrahedra < 1 || !out.neighborlist)
                    return;
                if (0) {
                    //DebugBdf(out, "tetgen.bdf", nullptr);
                }
                tetFlag.resize(out.numberoftetrahedra, 0);
                std::stack<int> elemStack;
                int nMarked = 0;
                for (int i = 0; i < out.numberoftetrahedra; i++) {
                    int index = i * 4;
                    for (int j = 0; j < 4; j++) {
                        int eAdj = out.neighborlist[index + j];
                        //outer tets.
                        if (eAdj < 0) {
                            tetFlag[i] = 1;
                            ++nMarked;
                            elemStack.push(i);
                            break;
                        }
                    }
                }
                //whether tet contains boundaryfaces.
                DInt3Array boundaryFace(in.numberoffacets);
                for (size_t i = 0; i < in.numberoffacets; i++) {
                    const auto& p = in.facetlist[i].polygonlist[0];
                    boundaryFace[i] = Int3(p.vertexlist[0], p.vertexlist[1], p.vertexlist[2]);
                }

                Int3HArray tetFacesProxy(boundaryFace);

                while (nMarked < out.numberoftetrahedra) {
                    int nNewMarked = 0;
                    DIntArray elemStackNext;
                    while (!elemStack.empty()) {
                        int ei = elemStack.top();
                        assert(tetFlag[ei] != 0);
                        int index = ei * 4;
                        elemStack.pop();
                        auto tet = out.tetrahedronlist + ei * 4;
                        for (int i = 0; i < 4; i++) {
                            int eAdj = out.neighborlist[index + i];
                            if (eAdj < 0 || tetFlag[eAdj] != 0)
                                continue;
                            //cross boundary face. face-topo may change!
                            int pos = tetFacesProxy.find(Int3(tet[TET_FACE_VERTEX[i][0]], tet[TET_FACE_VERTEX[i][1]], tet[TET_FACE_VERTEX[i][2]]));
                            if (pos < 0) {
                                elemStack.push(eAdj);
                                tetFlag[eAdj] = tetFlag[ei];
                                nMarked++;
                                ++nNewMarked;
                            }
                            else {
                                elemStackNext.push_back(eAdj);
                                tetFlag[eAdj] = -tetFlag[ei];
                                nMarked++;
                                ++nNewMarked;
                            }
                        }
                    }
                    if (nNewMarked < 1)
                        break;
                    Unique(elemStackNext);
                    for (auto it : elemStackNext) {
                        elemStack.push(it);
                    }
                }
            };

            std::stringstream os;
            ScopedRedirect redirect;
            tetgenio::gPrintFun = MyPrintFun;
            tetgenio tetIn, out, addin;
            Json::Value jsonObj;
            tetgenbehavior b;
            {
                if (!SetByFile(jsonObj, finput))
                    return -1;

                std::string meshCtrl, meshdata, addmeshdata;

                if (!GetValue(jsonObj, "tetgenInputMesh", meshdata))
                {
                    return -1;
                }

                if (HasKey(jsonObj, "behavior"))
                {
                    meshCtrl = Convert2TetGenCommand(jsonObj["behavior"], tetIn.numberoftetrahedra);
                }
                else if (!GetValue(jsonObj, "behavior4Tetgen", meshCtrl))
                {
                    return false;
                }
                b.parse_commandline((char*)(meshCtrl.c_str()));

                if (!tetgeniobyfile(tetIn, meshdata))
                {
                    if (!tetgeniobystring(tetIn, meshdata))
                    {
                        return false;
                    }
                }

                gRefineData.mDfMax = -1;

                double dfmax;
                if (tGetValue(jsonObj, "meshSizeMax", dfmax) && dfmax > 0.0) {
                    gRefineData.mDfMax= dfmax;
                    tetIn.tetunsuitable = IsTetraNeedRefine;
                }

                std::string refinementsFile;
                if (GetValue(jsonObj, "refinements", refinementsFile))
                {
                    auto lmdReadRefine = [](std::istream& mIO)
                    {
                        size_t num;
                        mIO.read((char*)(&num), sizeof(size_t));
                        gRefineData.mValue.resize(num);
                        mIO.read((char*)gRefineData.mValue.data(), sizeof(float)*gRefineData.mValue.size());

                        mIO.read((char*)gRefineData.sizeData, 4 * sizeof(double));
                        mIO.read((char*)gRefineData.mSize, 3 * sizeof(int32_t));
                        mIO.read((char*)gRefineData.mStart, 3 * sizeof(double));
                        mIO.read((char*)gRefineData.mStep, 3 * sizeof(double));
                        mIO.read((char*)gRefineData.mStepInv, 3 * sizeof(double));
                    };

                    std::ifstream mIO(refinementsFile.c_str(), std::ios_base::in | std::ios_base::binary);
                    if (!mIO) {
                        std::string  strOut;

                        if (Base64::Decode(refinementsFile, &strOut))
                        {
                            std::stringstream mIO;
                            mIO.write(strOut.data(), strOut.size());
                            lmdReadRefine(mIO);
                            tetIn.tetunsuitable = IsTetraNeedRefine;
                        }
                    }
                    else
                    {
                        lmdReadRefine(mIO);
                        tetIn.tetunsuitable = IsTetraNeedRefine;
                    }
                }

                if (GetValue(jsonObj, "addInputMesh", addmeshdata))
                {
                    tetgeniobystring(addin, addmeshdata);
                }

            }
            tetrahedralize(&b, &tetIn, &out, &addin, nullptr);
            if (out.numberofpoints < 4 || out.numberoftetrahedra < 1)
                return -1;

            int markTetra = 1;
            tGetValue(jsonObj, "markTetra", markTetra);

            
            std::vector<int> tetFlag;
            if (markTetra > 0) {
                MarkTetra(tetIn, out, tetFlag);
            }

            int pos = fout.find_last_of('.');
            if (fout.substr(pos) == ".json") {
                std::string strData;
                if (!tetgenio2string(out, strData, markTetra > 0? &tetFlag:nullptr))
                {
                    return false;
                }

                Json::Value result;
                result["meshResult"] = strData;
                std::ofstream out(fout);
                Json::FastWriter wr;
                out << wr.write(result);

            }
            else if (!tetgenio2file(out, fout, markTetra > 0? &tetFlag:nullptr))
                return false;
        }
        catch (...) {
            return -1;
        }

        return 0;
    }

    int MeshGen(const std::string& finput, const std::string& fout)
    {
        typedef ArrayHashProxy<DInt3Array, HashInt3> Int3HArray;

        try {
            Json::Value jsonObj, jsonOut;
            if (!SetByFile(jsonObj, finput))
                return -1;

            tetgenio tetIn, out, inadd;

            auto ImportMesh = [&]()
            {
                std::vector<double> points;
                if (!GetData(jsonObj, "meshInput.points", points))
                {
                    return false;
                }
                tetIn.numberofpoints = points.size() / 3;
                tetIn.pointlist = new REAL[tetIn.numberofpoints * 3];
                memcpy(tetIn.pointlist, points.data(), points.size() * sizeof(double));

                if (GetData(jsonObj, "refinements.points", points))
                {
                    // bgmin.numberofpoints = points.size() / 4;
                    // bgmin.pointlist = new REAL[bgmin.numberofpoints * 3];
                    // bgmin.numberofpointmtrs = 1;
                    // bgmin.pointmtrlist = new REAL[bgmin.numberofpoints * bgmin.numberofpointmtrs];
                    //
                    //
                    // for (int i = 0; i < bgmin.numberofpoints; i++)
                    // {
                    //   for (int j = 0; j < 3; j++)
                    //   {
                    //     tetIn.pointlist[i * 3 + j] = points[i * 4 + j];
                    //   }
                    //   bgmin.pointmtrlist[i] = points[i * 4 + 3];
                    // }
                }

                {
                    std::vector<int>    triangles, quad;
                    GetData(jsonObj, "meshInput.triangles", triangles);
                    GetData(jsonObj, "meshInput.quadrilaterals", quad);
                    tetIn.numberoffacets = triangles.size() / 4 + quad.size() / 5;
                    if (tetIn.numberoffacets > 0) {
                        tetIn.facetlist = new tetgenio::facet[tetIn.numberoffacets];
                        tetIn.facetmarkerlist = new int[tetIn.numberoffacets];
                        int fi = 0;
                        auto ntri = triangles.size() / 4;
                        for (size_t i = 0; i < ntri; i++)
                        {
                            tetgenio::facet *f = &tetIn.facetlist[fi];
                            tetgenio::init(f);
                            f->numberofholes = 0;
                            f->numberofpolygons = 1;
                            f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
                            tetgenio::polygon *p = &f->polygonlist[0];
                            tetgenio::init(p);
                            p->numberofvertices = 3;
                            p->vertexlist = new int[p->numberofvertices];
                            for (int j = 0; j < p->numberofvertices; j++)
                            {
                                p->vertexlist[j] = triangles[i * 4 + j];
                            }

                            tetIn.facetmarkerlist[fi] = 203;
                            ++fi;
                        }

                        int nquad = quad.size() / 5;
                        for (size_t i = 0; i < nquad; i++)
                        {
                            tetgenio::facet *f = &tetIn.facetlist[fi];
                            tetgenio::init(f);
                            f->numberofholes = 0;
                            f->numberofpolygons = 1;
                            f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
                            tetgenio::polygon *p = &f->polygonlist[0];
                            tetgenio::init(p);
                            p->numberofvertices = 4;
                            p->vertexlist = new int[p->numberofvertices];
                            for (int j = 0; j < p->numberofvertices; j++)
                            {
                                p->vertexlist[j] = quad[i * 5 + j];
                            }

                            tetIn.facetmarkerlist[fi] = 204;
                            ++fi;
                        }
                    }
                }

                {
                    std::vector<int> tets;
                    GetData(jsonObj, "meshInput.tetrahedrons", tets);
                    int ntet = tets.size() / 5;
                    tetIn.numberoftetrahedra = ntet;
                    if (ntet > 0) {
                        tetIn.tetrahedronlist = new int[tetIn.numberoftetrahedra * 4];
                        if (!tetIn.tetrahedronlist) {
                            return false;
                        }
                        int index = 0;
                        for (size_t i = 0; i < ntet; i++) {
                            for (int j = 0; j < 4; j++) {
                                tetIn.tetrahedronlist[index++] = tets[i * 5 + j];
                            }
                        }
                    }
                }

                return true;
            };

            auto  MarkTetra = [](const tetgenio& in, const tetgenio& out, DIntArray& tetFlag) {
                if (out.numberoftetrahedra < 1 || !out.neighborlist)
                    return;
                if (0) {
                    //DebugBdf(out, "tetgen.bdf", nullptr);
                }
                tetFlag.resize(out.numberoftetrahedra, 0);
                std::stack<int> elemStack;
                int nMarked = 0;
                for (int i = 0; i < out.numberoftetrahedra; i++) {
                    int index = i * 4;
                    for (int j = 0; j < 4; j++) {
                        int eAdj = out.neighborlist[index + j];
                        //outer tets.
                        if (eAdj < 0) {
                            tetFlag[i] = 1;
                            ++nMarked;
                            elemStack.push(i);
                            break;
                        }
                    }
                }
                //whether tet contains boundaryfaces.
                DInt3Array boundaryFace(in.numberoffacets);
                for (size_t i = 0; i < in.numberoffacets; i++) {
                    const auto& p = in.facetlist[i].polygonlist[0];
                    boundaryFace[i] = Int3(p.vertexlist[0], p.vertexlist[1], p.vertexlist[2]);
                }

                Int3HArray tetFacesProxy(boundaryFace);

                while (nMarked < out.numberoftetrahedra) {
                    int nNewMarked = 0;
                    DIntArray elemStackNext;
                    while (!elemStack.empty()) {
                        int ei = elemStack.top();
                        assert(tetFlag[ei] != 0);
                        int index = ei * 4;
                        elemStack.pop();
                        auto tet = out.tetrahedronlist + ei * 4;
                        for (int i = 0; i < 4; i++) {
                            int eAdj = out.neighborlist[index + i];
                            if (eAdj < 0 || tetFlag[eAdj] != 0)
                                continue;
                            //cross boundary face. face-topo may change!
                            int pos = tetFacesProxy.find(Int3(tet[TET_FACE_VERTEX[i][0]], tet[TET_FACE_VERTEX[i][1]], tet[TET_FACE_VERTEX[i][2]]));
                            if (pos < 0) {
                                elemStack.push(eAdj);
                                tetFlag[eAdj] = tetFlag[ei];
                                nMarked++;
                                ++nNewMarked;
                            }
                            else {
                                elemStackNext.push_back(eAdj);
                                tetFlag[eAdj] = -tetFlag[ei];
                                nMarked++;
                                ++nNewMarked;
                            }
                        }
                    }
                    if (nNewMarked < 1)
                        break;
                    Unique(elemStackNext);
                    for (auto it : elemStackNext) {
                        elemStack.push(it);
                    }
                }
            };

            auto  ExportMesh = [&]()
            {
                std::vector<int> tetFlag;
                MarkTetra(tetIn, out, tetFlag);

                std::vector<double> points;

                points.resize(out.numberofpoints * 3);
                memcpy(points.data(), out.pointlist, out.numberofpoints * 3 * sizeof(double));

                if (!WriteData(jsonOut, "meshOutput.points", points))
                {
                    return false;
                }

                std::vector<int>    tets;
                int Ntet = 0;
                for (size_t i = 0; i < out.numberoftetrahedra; i++)
                {
                    if (tetFlag[i] < 0)
                        continue;
                    auto tet = out.tetrahedronlist + i * 4;
                    for (int j = 0; j < 4; j++) {
                        tets.push_back(tet[j]);
                    }
                    tets.push_back(304);
                    Ntet++;
                }
                if (!WriteData(jsonOut, "meshOutput.tetrahedrons", tets))
                {
                    return false;
                }

                return true;
            };

            if (!ImportMesh())
                return -1;
            if (0) {
                DebugBdf(tetIn, "tetgenIn.bdf", nullptr);
            }
            auto meshCtrl = Convert2TetGenCommand(jsonObj, tetIn.numberoftetrahedra);
            tetgenbehavior b;
            b.parse_commandline((char*)meshCtrl.c_str());
            tetrahedralize(&b, &tetIn, &out, inadd.numberofpoints > 0 ? &inadd : nullptr);
            if (out.numberofpoints < 4 || out.numberoftetrahedra < 1)
                return -1;
            if (0) {
                DebugBdf(out, "tetgenOut.bdf", nullptr);
            }

            // ouput
            if (!ExportMesh())
                return -1;

            if (!SetToFile(jsonOut, fout))
                return false;

        }
        catch (...) {
            return -1;
        }
        return 0;
    }
}