#include "TemplateMesher.h"
#include "misc/miscFunctions.hpp"
#include <TopTools_ListOfShape.hxx>
#include <TopTools_IndexedMapOfShape.hxx>
#include <Bnd_Box.hxx>
#include <BRepBndLib.hxx>
#include <TopExp_Explorer.hxx>
#include <GProp_GProps.hxx>
#include <BRepGProp.hxx>
#include <fstream>

#include "topology/ShapeAlgorithms.h"
#include "mesh/GmshMesh.hpp"
#include "mesh/MeshFactory.hpp"
#include "mesh/MeshToVtkWriter.hpp"
#include "misc/Timer.hpp"

namespace nglib {
#include "nglib.h"
}

TemplateMesher::TemplateMesher(TopoDS_Shape c)
{
    meshDensityFactor = 1.0;
    compound = c;

    //calculate minimal and maximal distances
    Bnd_Box box;
    box.SetGap(0.0);
    BRepBndLib::Add(compound, box);

    Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
    box.Get(xmin, ymin, zmin, xmax, ymax, zmax);

    double maxdist = std::sqrt( pow( (xmax - xmin),2) + pow((ymax - ymin), 2) + pow((zmax - zmin), 2));

    double mindist = maxdist;

    for (TopExp_Explorer edgeExplorer(compound, TopAbs_ShapeEnum::TopAbs_FACE); edgeExplorer.More(); edgeExplorer.Next())
    {
        TopoDS_Shape thisEdge = edgeExplorer.Current();
        GProp_GProps massProps;
        BRepGProp::LinearProperties(thisEdge, massProps);
        if (massProps.Mass() < mindist)
            mindist = massProps.Mass();
    }

    globalMinH = std::max(mindist, maxdist/500);
    globalMaxH = maxdist/10;
    meshGrading = 0.4;
}

TemplateMesher::~TemplateMesher()
{
    delete meshPtr;
}

void TemplateMesher::writeMesh()
{
    using namespace nglib;
    Ng_Mesh* occ_mesh = (Ng_Mesh*)meshPtr;
    // writing mesh for cenos
    misc::createFolder(projectPath + "\\geometry");
    std::ofstream outfile(projectPath + "\\geometry\\meshFile.msh");
    outfile.precision(6);

    outfile << "$MeshFormat\n";
    outfile << (float)2.0 << " "
        << (int)0 << " "
        << (int)sizeof(double) << "\n";
    outfile << "$EndMeshFormat\n";


    outfile << "$PhysicalNames\n";
    outfile << emap.size() << "\n";


    for (auto& ent : emap.getEntities())
        outfile << ent->getId() << " " << ent->getName() << "\n";

    outfile << "$EndPhysicalNames\n";

    outfile << "$Nodes\n";
    int NrPoints = Ng_GetNP(occ_mesh);
    outfile << NrPoints << "\n";
    for (int j = 1; j <= NrPoints; j++)
    {
        double p[3];
        Ng_GetPoint(occ_mesh, j, p);
        outfile << j << " " << p[0] << " " << p[1] << " " << p[2] << "\n";
    }
    outfile << "$EndNodes\n";

    outfile << "$Elements\n";
    int NrSurfEl = Ng_GetNSE(occ_mesh);
    int NrVolEl = Ng_GetNE(occ_mesh);
    outfile << NrSurfEl + NrVolEl << "\n";
    for (int j = 1; j <= NrSurfEl; j++)
    {
        int nodes[3];
        Ng_GetSurfaceElement(occ_mesh, j, nodes);
        // index / element type (trig = 2) / nr of tags (2) / physicl tag / elementary tag / nodes...
        int physicalTag = Cenos_GetSurfaceElementIndex(occ_mesh, j);
        int entityId = emap.getEntityByFaceTopologyId(physicalTag)->getId();
        outfile << j << " " << "2" << " " << "2" << " " << entityId << " " << entityId << " ";
        outfile << nodes[0] << " " << nodes[1] << " " << nodes[2] << "\n";
    }

    for (int j = 1; j <= Ng_GetNE(occ_mesh); j++)
    {
        int nodes[4];
        Ng_GetVolumeElement(occ_mesh, j, nodes);
        // index / element type (tetra = 4) / nr of tags (2) / physicl tag / elementary tag / nodes...
        int physicalTag = Cenos_GetVolumeElementIndex(occ_mesh, j);
        int entityId = emap.getEntityBySolidTopologyId(physicalTag)->getId();
        outfile << j + NrSurfEl << " " << "4" << " " << "2" << " " << entityId << " " << entityId << " ";
        outfile << nodes[0] << " " << nodes[1] << " " << nodes[2] << " " << nodes[3] << "\n";
    }
    outfile << "$EndElements\n";

    outfile.close();
}


void TemplateMesher::calculateEntityMap()
{
    using namespace nglib;

    TopoDS_Shape* final_shape = &compound;
    Ng_OCC_Geometry* occ_geom = Cenos_OCC_ShapeToGeometry((Ng_TopoDS_Shape*)final_shape);
    TopTools_IndexedMapOfShape FMap;
    Ng_OCC_TopTools_IndexedMapOfShape* occ_fmap = (Ng_OCC_TopTools_IndexedMapOfShape*)&FMap;

    Ng_Result ng_res;

    ng_res = Ng_OCC_GetFMap(occ_geom, occ_fmap);

    if (!FMap.Extent() or (ng_res !=Ng_Result::NG_OK))
    {
        std::cout << "Error retrieving Face map...." << std::endl;
        return;
    }

    int entityIndex = 0;
    for (int i = 1; i <= FMap.Extent(); i++)
    {
        TopoDS_Shape OCCface;
        OCCface = FMap.FindKey(i);
        for (auto& face : faceMap)
        {
            if (ShapeAlgorithms::have_common(OCCface, face.second))
            {
                //check if this face index is already created
                //if it is created, use it for faceIdsMap with current i
                std::shared_ptr<Entity> ent = emap.getEntityByName(face.first);
                if (ent  != nullptr)
                {
                    ent->addTopologyId(i);
                    break;
                }
                else
                {
                    entityIndex++;
                    Entity entNew(entityIndex, face.first);
                    entNew.addTopologyId(i);
                    entNew.setEntityType(Entity::AREA);
                    emap.addEntity(entNew);
                }
            }
        }

        std::shared_ptr<Entity> ent = emap.getEntityByFaceTopologyId(i);
        if (ent == nullptr)
        {
            entityIndex++;
            Entity entNew(entityIndex, "_Face" + std::to_string(i));
            entNew.addTopologyId(i);
            entNew.setEntityType(Entity::AREA);
            emap.addEntity(entNew);
        }
    }

    //face map for recognizing shapes
    TopTools_IndexedMapOfShape SoMap;
    Ng_OCC_TopTools_IndexedMapOfShape* occ_somap = (Ng_OCC_TopTools_IndexedMapOfShape*)&SoMap;

    ng_res = Cenos_OCC_GetSoMap(occ_geom, occ_somap);

    if (!SoMap.Extent() or (ng_res != Ng_Result::NG_OK))
    {
        std::cout << "Error retrieving Solid map...." << std::endl;
        return;
    }

    for (int i = 1; i <= SoMap.Extent(); i++)
    {

        TopoDS_Shape OCCsolid;
        OCCsolid = SoMap.FindKey(i);
        for (auto& solid : solidMap)
        {
            // skip air, as first all other solids have to be established
            if (solid.first == "air")
                continue;

            if (ShapeAlgorithms::have_common(OCCsolid, solid.second))
            {
                std::shared_ptr<Entity> ent = emap.getEntityByName(solid.first);
                if (ent != nullptr)
                {
                    ent->addTopologyId(i);
                    break;
                }
                else
                {
                    entityIndex++;
                    Entity entNew(entityIndex, solid.first);
                    entNew.addTopologyId(i);
                    entNew.setEntityType(Entity::VOLUME);
                    emap.addEntity(entNew);
                }
            }
        }

        //now add air, if this OCCSolid does not belong to any shape
        std::shared_ptr<Entity> ent = emap.getEntityBySolidTopologyId(i);
        if (ent == nullptr)
        {
            ent = emap.getEntityByName("air");
            if (ent != nullptr)
            {
                ent->addTopologyId(i);
            }
            else
            {
                entityIndex++;
                Entity entNew(entityIndex, "air");
                entNew.addTopologyId(i);
                entNew.setEntityType(Entity::VOLUME);
                emap.addEntity(entNew);
            }
        }
    }

}

void TemplateMesher::refineSubMeshes()
{
    for (auto sm : edgeSubMeshes)
    {
        for (TopoDS_Shape& edge : sm.shapeList)
        {
            using namespace nglib;
            Bnd_Box B;
            double p1[3], p2[3];
            BRepBndLib::Add(edge, B);
            B.Get(p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]);
            Ng_RestrictMeshSizeBox((Ng_Mesh*)meshPtr, p1, p2, globalMinH * sm.meshSizeFactor * meshDensityFactor);
        }
    }
}


void TemplateMesher::setMinMeshSize(double s)
{
    this->globalMinH = s;
}

void TemplateMesher::setMaxMeshSize(double s)
{
    this->globalMaxH = s;
}


void TemplateMesher::setMeshGrading(double s)
{
    if ( (s > 1.0) or (s < 0.0) )
    {
        this->meshGrading = 0.4;
    }
    else
        this->meshGrading = s;
}

void TemplateMesher::setRefinementOnEdges(TopTools_ListOfShape edgeList, double factor)
{
    edgeSubMeshes.push_back( SubMesh(edgeList, factor) );
}

json TemplateMesher::getMeshEntities()
{
    return meshEntities;
}


void TemplateMesher::setProjectPath(std::string p)
{
    projectPath = p;
}

void TemplateMesher::setFaceMap(std::map<std::string, TopTools_ListOfShape> fm)
{
    faceMap = fm;
}

void TemplateMesher::setSolidMap(std::map<std::string, TopTools_ListOfShape> sm)
{
    solidMap = sm;
}

void TemplateMesher::generateMesh()
{
    using namespace nglib;

    TopoDS_Shape* final_shape = &compound;
    Ng_OCC_Geometry* occ_geom = Cenos_OCC_ShapeToGeometry((Ng_TopoDS_Shape*)final_shape);
    Ng_Mesh* occ_mesh;
    Ng_Meshing_Parameters mp;

    // Result of Netgen Operations
    Ng_Result ng_res;

    // Initialise the Netgen Core library
    Ng_Init();
    
    //Redirect after init!!! Ng_Init sets cout to std::cout
    std::ofstream fileout(projectPath + "\\mesh_log.log");
    if (fileout.is_open())
    {
        std::ofstream* ptr_fileout = &fileout;
        Cenos_RedirectCout(ptr_fileout);
    }
    else
    {
        std::cout << "WARNING: could not redirect meshing output to log file, will be written to console" << endl;;
    }

    occ_mesh = Ng_NewMesh();
    meshPtr = occ_mesh;

    mp.uselocalh = 1;
    mp.elementsperedge = 2.0;
    mp.elementspercurve = 2.0;
    mp.maxh = globalMaxH * meshDensityFactor;
    mp.minh = globalMinH * meshDensityFactor;
    mp.grading = meshGrading;
    mp.closeedgeenable = 0;
    mp.closeedgefact = 1.0;
    mp.optsurfmeshenable = 1;

    std::cout << "Setting Local Mesh size....." << std::endl;
    Ng_OCC_SetLocalMeshSize(occ_geom, occ_mesh, &mp);
    std::cout << "Local Mesh size successfully set....." << std::endl;

    //refineSubMeshes();

    std::cout << "Creating Edge Mesh....." << std::endl;
    ng_res = Ng_OCC_GenerateEdgeMesh(occ_geom, occ_mesh, &mp);


    if (ng_res != NG_OK)
    {
        Ng_DeleteMesh(occ_mesh);
        std::cout << "Error creating Edge Mesh.... Aborting!!" << std::endl;
        return;
    }
    else
    {
        std::cout << "Edge Mesh successfully created....." << std::endl;
        std::cout << "Number of points = " << Ng_GetNP(occ_mesh) << std::endl;
    }

    std::cout << "Creating Surface Mesh....." << std::endl;

    try
    {
        ng_res = Ng_OCC_GenerateSurfaceMesh(occ_geom, occ_mesh, &mp);
    }
    catch (std::exception& e)
    {
        std::cout << "Exception during surface meshing " << e.what() << std::endl;
        return;
    }

    if (ng_res != NG_OK)
    {
        Ng_DeleteMesh(occ_mesh);
        std::cout << "Error creating Surface Mesh..... Aborting!!" << std::endl;
        return;
    }
    else
    {
        std::cout << "Surface Mesh successfully created....." << std::endl;
        std::cout << "Number of points = " << Ng_GetNP(occ_mesh) << std::endl;
        std::cout << "Number of surface elements = " << Ng_GetNSE(occ_mesh) << std::endl;
    }

    std::cout << "Creating Volume Mesh....." << std::endl;

    try
    {
        ng_res = Ng_GenerateVolumeMesh(occ_mesh, &mp);
    }
    catch (std::exception& e)
    {
        std::cout << "Exception during volume meshing " << e.what() << std::endl;
        return;
    }

    if (ng_res != NG_OK)
    {
        Ng_DeleteMesh(occ_mesh);
        std::cout << "Error creating Volume Mesh..... Aborting!!" << std::endl;
        return;
    }
    else
    {
        std::cout << "Volume Mesh successfully created....." << std::endl;
        std::cout << "Number of points = " << Ng_GetNP(occ_mesh) << std::endl;
        std::cout << "Number of volume elements = " << Ng_GetNE(occ_mesh) << std::endl;
    }

    misc::createFolder(projectPath + "\\geometry");

    calculateEntityMap();

    //writing via Mesh class is preferred, as it provides access to getMesh Entities
    //writeMesh();

    std::shared_ptr<Mesh> mesh = std::make_shared<GmshMesh>(occ_mesh, emap);

    mesh->processMesh();

    std::ofstream meshOutStreamFile;

    meshOutStreamFile.open(projectPath + "\\geometry\\meshFile.msh");
    mesh->writeMesh(meshOutStreamFile);
    meshOutStreamFile.close();

    fileout.close();
}

void TemplateMesher::setMeshDensity(std::string md)
{
    if (md == "rough")
        meshDensityFactor = 1.5;
    else if (md == "fine")
        meshDensityFactor = 0.55;
    else
        meshDensityFactor = 1.0;
}

TemplateMesher::SubMesh::SubMesh(TopTools_ListOfShape sl, double ms): shapeList(sl), meshSizeFactor(ms)
{
}