#include "MeshGenerator.h"
#include <filesystem>
#include <stdexcept>
#include <Bnd_Box.hxx>
#include <BRepBndLib.hxx>
#include <TopExp_Explorer.hxx>
#include <GProp_GProps.hxx>
#include <BRepGProp.hxx>
#include <BRepAlgoAPI_Fuse.hxx>
#include <BRepBuilderAPI_Transform.hxx>
#include <BRepAlgoAPI_Cut.hxx>
#include <BRepAlgoAPI_BuilderAlgo.hxx>
#include <BRep_Tool.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Vertex.hxx>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <Standard_Type.hxx>
#include <Geom_Line.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <BRepTools.hxx>
#include "topology/ShapeAlgorithms.h"
#include "topology/ShapePartition.h"
#include "misc/miscFunctions.hpp"
#include "mesh/MeshEntity.h"
#include "cenos_exception.h"
#include <TopExp.hxx>
#include "NetgenWrapper/NetgenWrapper.h"


MeshGenerator::~MeshGenerator()
{
}

MeshGenerator::MeshGenerator(TopoDS_Shape part): mesh_parameters(std::make_shared<MeshParameters>(part)), stopper(nullptr)
{
    //partition will be used for meshing
    partition = part;
  
    ShapeType shape_type = TopoEntity::getShapeType(part);
    if (shape_type == ShapeType::SOLID)
    {
        mesh_mode = M_3D;
        mesh_shape_type = ShapeType::SOLID;
    }
    else if (shape_type == ShapeType::FACE)
    {
        mesh_mode = M_2D;
        mesh_shape_type = ShapeType::FACE;
    }
    else if (shape_type == ShapeType::EDGE)
    {
        mesh_mode = M_1D;
        mesh_shape_type = ShapeType::EDGE;
    }
    else
    {
        std::string message = "Cannot initialize Mesh Generator, cannot determine type of shape";
        throw cenos_exception(message);
    }
    main_ng_mesh = std::make_shared< NetgenMesh>();
}


MeshGenerator::MeshGenerator(GeometryData_ geom_data): stopper(nullptr)
{
    partition = geom_data->getShape();
    mesh_parameters = std::make_shared<MeshParameters>(partition);

    for (auto ent : geom_data->getEntities())
    {
        if (ent->isEnabled())
            addTopoEntity(ent, true); // force adding entity, this saves time. geometryData should have non-overlapping entities
    }

    relations = geom_data->getRelations();

    ShapeType shape_type = TopoEntity::getShapeType(partition);
    if (shape_type == ShapeType::SOLID)
    {
        mesh_mode = M_3D;
        mesh_shape_type = ShapeType::SOLID;
    }
    else if (shape_type == ShapeType::FACE)
    {
        mesh_mode = M_2D;
        mesh_shape_type = ShapeType::FACE;
    }
    else if (shape_type == ShapeType::EDGE)
    {
        mesh_mode = M_1D;
        mesh_shape_type = ShapeType::EDGE;
    }
    else
    {
        std::string message = "Cannot initialize Mesh Generator, cannot determine type of shape";
        throw cenos_exception(message);
    }
    main_ng_mesh = std::make_shared< NetgenMesh>();
}

MeshGenerator::MeshGenerator(const MeshGenerator& that)
{
    wDir = that.wDir;
    partition = that.partition;
    meshing_data = that.meshing_data;
    mesh_parameters = that.mesh_parameters;
    stopper = that.stopper;
}


void MeshGenerator::setWdir(std::string w)
{
    wDir = w;
}

void MeshGenerator::setProgressIndicator(std::shared_ptr<ProgressIndicator> pi)
{
    progress = pi;
}


//returns true if succesfully added, false if such entity exists
bool MeshGenerator::addTopoEntity(std::shared_ptr<TopoEntity> te, bool force_add)
{
    if (ShapeAlgorithms::getHighestShapeType(te->getShape()) == TopAbs_EDGE)
    {
        TopoDS_Edge edge = TopoDS::Edge(te->getShape());
        if (BRep_Tool::Degenerated(edge))
        {
            // Ignoring degenerate TopoEntity 
            return false;
        }
    }

    // if force_add, do not check for overlapping! this saves a lot of time 
    if (!force_add)
    {
        for (auto ent : entities)
        {
        if ((ShapeAlgorithms::have_equal_shapes(te->getShape(), ent->getShape()))
                && (ShapeAlgorithms::sameType(te->getShape(), ent->getShape())))
            {
                return false;
            }
        }
    }

    if (te->getId() == 0)
        te->setId(entities.size() + 1);
    entities.push_back(te);
    std::shared_ptr<MeshEntity> me = std::make_shared<MeshEntity>(te->getShape(), te->getName(), te->getId());
    me->setMeshParameters(mesh_parameters);
    meshing_data.addMeshEntity(me);

    return true;
}

void MeshGenerator::setParameters(std::shared_ptr<MeshParameters> mp)
{
    mesh_parameters = mp;
    for (auto ent : entities)
    { 
        meshing_data.setParameters(ent, mesh_parameters);
    }
}

/// <summary>
/// set parameters mp to entity and all subentities!
/// </summary>
/// <param name="mp"></param>
/// <param name="ent"></param>
void MeshGenerator::setParameters(std::shared_ptr<MeshParameters> mp, std::shared_ptr<TopoEntity> ent0)
{
    meshing_data.setParameters(ent0, mp);
    // set mp for faces (3D mesh) or edges (2D mesh)
    for (auto ent : entities)
    {
        if ( (ent->getDimension() == ent0->getDimension() - 1) and ent->isSubShapeOf(ent0->getShape()))
            meshing_data.setParameters(ent, mp);
    }

    // set mp for edges (3D mesh)
    for (auto ent : entities)
    {
        if ( (ent->getDimension() == ent0->getDimension() - 2) and ent->isSubShapeOf(ent0->getShape()))
            meshing_data.setParameters(ent, mp);
    }
}

void MeshGenerator::setTransform(std::shared_ptr<TopoEntity> ent, gp_Trsf transform)
{
    entity_transform[ent] = transform;
}

void MeshGenerator::setStopper(bool* b)
{
    stopper = b;
}

std::shared_ptr<Mesh> MeshGenerator::generateMesh()
{
    validateInput();

    // initialize netgen (set logging file)
    // must be done only once on starting of meshing
    InitializeNetgen();

    copyTransformedMeshData();

    NetgenGeometryWrapper ng_geom(&partition);

    auto mp = mesh_parameters;

    main_ng_mesh->SetMeshSizing(&ng_geom, mp);

    //TODO make meshing order!
    switch (mesh_mode) 
    {
        case Mode::M_1D:
        {
            for (auto ent : entities)
            {
                if (ent->getType() == ShapeType::EDGE)
                {
                    NetgenMesh_ edge_mesh = generate1DMesh(ent);
                    main_ng_mesh->MergeMesh(edge_mesh, ent->getId());
                }
            }
            break;
        }
        case Mode::M_2D:
        {
            if (progress != nullptr)
            {
                progress->setProgress(1);
            }
            
            int nr_of_faces = 0;
            for (auto ent : entities)
            {
                if (ent->getType() == ShapeType::FACE)
                    nr_of_faces++;
            }

            int i = 0;
            for (auto ent : entities)
            {
                if (ent->getType() == ShapeType::FACE)
                {
                    i++;
                    if (progress != nullptr)
                        progress->sendLog("Meshing face " + std::to_string(i) + " / " + std::to_string(nr_of_faces));

                    NetgenMesh_ face_mesh = generate2DMesh(ent);

                    main_ng_mesh->MergeMesh(face_mesh, ent->getId());
                    if (progress != nullptr)
                        progress->setProgress(100 * i / nr_of_faces);
                }
            }

            //copy missing edges
            for (auto ent : entities)
            {
                if (ent->getType() == ShapeType::EDGE)
                {
                    if (ent->getShape().Orientation() == TopAbs_Orientation::TopAbs_REVERSED)
                        meshing_data.getMesh(ent)->ReverseSegments();
                    main_ng_mesh->MergeMesh(meshing_data.getMesh(ent), ent->getId());
                }
            }

            break;
        }
        case Mode::M_3D:
        {
            //solid map for recognizing shapes

            TopTools_IndexedMapOfShape SolidMap = ng_geom.GetSolidMap();

            if (progress != nullptr)
            {
                progress->setProgress(1);
            }

            for (int i = 1; i <= SolidMap.Extent(); i++)
            {
                if (progress != nullptr)
                    progress->sendLog("Meshing solid " + std::to_string(i) + " / " + std::to_string(SolidMap.Extent()));
                bool match_found = false;
                TopoDS_Shape OCCsolid;
                OCCsolid = SolidMap.FindKey(i);

                for (auto s : entities)
                {
                    if (s->getType() == ShapeType::SOLID)
                    {
                        if (ShapeAlgorithms::have_equal_shapes(OCCsolid, s->getShape()))
                        {

                            solid_id[i] = s->getId();

                            NetgenMesh_ solid_mesh = generate3DMesh(s);
                            main_ng_mesh->MergeMesh(solid_mesh, i);
                            match_found = true;
                            break;
                        }
                    }
                }

                if (!match_found) {
                    throw cenos_exception("No matching entity found for OCCsolid_" + std::to_string(i));
                }

                if (progress != nullptr)
                    progress->setProgress(100 * i / (SolidMap.Extent()));
            }

            //copy missing edges

            //edge map for recognizing shapes
            TopTools_IndexedMapOfShape EdgeMap = ng_geom.GetEdgeMap();

            for (int i = 1; i <= EdgeMap.Extent(); i++)
            {
                bool match_found = false;
                TopoDS_Shape OCCedge;
                OCCedge = EdgeMap.FindKey(i);

                for (auto s : entities)
                {
                    if (s->getType() == ShapeType::EDGE)
                    {
                        if (ShapeAlgorithms::have_equal_shapes(OCCedge, s->getShape()))
                        {
                            edge_id[i] = s->getId();
                            if (s->getShape().Orientation() == TopAbs_Orientation::TopAbs_REVERSED)
                                meshing_data.getMesh(s)->ReverseSegments();
                            main_ng_mesh->MergeMesh(meshing_data.getMesh(s), i);
                            match_found = true;
                            break;
                        }
                    }
                }

                if (!match_found) {
                    if (!BRep_Tool::Degenerated(TopoDS::Edge(OCCedge))) {
                        // Degenerated edges might not have any mesh, so that scenario is ok.
                        throw cenos_exception("No matching entity found for OCCedge_" + std::to_string(i));
                    }
                }
            }

            //face map for recognizing shapes
            TopTools_IndexedMapOfShape FaceMap = ng_geom.GetFaceMap();

            for (int i = 1; i <= FaceMap.Extent(); i++)
            {
                bool match_found = false;
                TopoDS_Shape OCCface;
                OCCface = FaceMap.FindKey(i);
                for (auto s : entities)
                {
                    if (s->getType() == ShapeType::FACE)
                    {
                        if (ShapeAlgorithms::have_equal_shapes(OCCface, s->getShape()))
                        {
                            face_id[i] = s->getId();
                            if (s->getShape().Orientation() == TopAbs_Orientation::TopAbs_REVERSED)
                                meshing_data.getMesh(s)->ReverseSegments();
                            main_ng_mesh->MergeMesh(meshing_data.getMesh(s), i);
                            match_found = true;
                            break;
                        }
                    }
                }
                if (!match_found) {
                    throw cenos_exception("No matching entity found for OCCface_" + std::to_string(i));
                }
            }

            break;
        }
        default:
        {
            std::string message = "Could not determine mesh mode. Please check the passed geometry";
            throw cenos_exception(message);
        }
    }

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    std::shared_ptr<Mesh> generated_mesh =  copyMeshData(main_ng_mesh);

    FinalizeNetgen();

    return generated_mesh;
}

std::shared_ptr<Mesh> MeshGenerator::refineMesh()
{
    NetgenGeometryWrapper ng_geom(partition);

    if (mesh == nullptr)
    {
        std::string message = "Cannot refine mesh. Mesh to refine is not specified.";
        throw cenos_exception(message);
    }

    NetgenMesh_ main_mesh = mesh->getNetgenMesh();

    InitializeNetgen();

    int i = 0;
    for (auto el : mesh->getElements())
    {
        if (el->getDimension() != 3)
            continue;
        i++;
        if (el->getMark() == 1)
            main_mesh->SetElementRefinement(i, true);
        else
            main_mesh->SetElementRefinement(i, false);
    }


    main_mesh->Refine();

    std::shared_ptr<Mesh> refined_mesh = copyMeshData(main_mesh);

    FinalizeNetgen();

    return refined_mesh;
}


void MeshGenerator::setMesh(std::shared_ptr<Mesh> m )
{
    mesh = m;
}

//
//
//   Private methods
//

void writeDebugFiles(std::string wDir, TopoDS_Shape sh, std::string name) {
    misc::createFolder(wDir + "/debug");
    ShapeAlgorithms::writeBrepFile(sh, wDir + "/debug/" + name);
}

void writeDebugFiles(std::string wDir, NetgenMesh_ mesh, std::string name) {
    misc::createFolder(wDir + "/debug");
    std::string fn = wDir + "/debug/" + name;
    mesh->SaveMesh(fn.c_str());
}

std::shared_ptr<Mesh> MeshGenerator::copyMeshData(NetgenMesh_ ng_mesh)
{
    std::shared_ptr<Mesh> copied_mesh = std::make_shared<Mesh>();
    copied_mesh->setMeshScale(1);

    /*****************************/
    /*  ADDING PHYSICAL ENTITIES */
    /*****************************/
    for (auto ent : entities)
    {
        std::shared_ptr<MeshEntity> e0 = std::make_shared<MeshEntity>(ent->getShape(), ent->getName(), ent->getId());
        e0->setEnabled(ent->isEnabled());
        copied_mesh->addEntity(e0);
    }

    /*****************/
    /*  ADDING NODES */
    /*****************/
    int nNodes = ng_mesh->GetNP();
    for (int i = 1; i <= nNodes; i++)
    {
        double p[3];
        ng_mesh->GetPoint(i, p);
        std::shared_ptr<Node> n0 = std::make_shared<Node>(i, p[0], p[1], p[2]);
        copied_mesh->addNode(n0);
    }

    /********************/
    /*  ADDING ELEMENTS */
    /********************/
    int nEdgeEl = ng_mesh->GetNSeg();
    int nFaceEl = ng_mesh->GetNSE();
    int nVolEl = ng_mesh->GetNE();
    int nElems = nEdgeEl + nFaceEl + nVolEl;

    //edge elements
    for (int i = 1; i <= nEdgeEl; i++)
    {

        int en[6] = {}; //element nodes
        int matnr = 0;
        ng_mesh->GetSegment(i, en); 

        ElementTypes elType = ElementTypes::LINE;
        std::shared_ptr<Element> e0 = Element::createElement(elType);

        std::vector<std::shared_ptr<Node>> nodeList;
        for (int j = 0; j < e0->getTypeNodeCount(); j++)
        {
            nodeList.push_back(copied_mesh->getNode(en[j]));
        }

        e0->setNodes(nodeList);
        e0->setId(i);
        int physicalTag = ng_mesh->GetEdgeIndexOfSegment(i);
        std::shared_ptr<MeshEntity> ent;
        if (edge_id.find(physicalTag) != edge_id.end())
            ent = copied_mesh->getEntityById(edge_id[physicalTag]);
        else
            ent = copied_mesh->getEntityById(physicalTag);

        if (ent == nullptr)
        {
            std::string message = "Failed to copy mesh data. Please check the input geometry! (getEntityById returned nullptr)";
            LOG(ERROR) << physicalTag;
            throw cenos_exception(message);
        }

        e0->setEntity(ent);

        copied_mesh->addElement(e0);
    }

    //surface elements
    for (int i = 1; i <= nFaceEl; i++)
    {
        int en[6]; //element nodes

        ElementTypes elType = ng_mesh->GetSurfaceElement(i, en);
        if (elType != ElementTypes::TRI and elType != ElementTypes::QUAD)
            throw(cenos_exception("Only TRIG and QUAD elements supported for 2D mesh from Netgen!"));

        std::shared_ptr<Element> e0 = Element::createElement(elType);

        std::vector<std::shared_ptr<Node>> nodeList;
        for (int j = 0; j < e0->getTypeNodeCount(); j++)
        {
            nodeList.push_back(copied_mesh->getNode(en[j]));
        }

        e0->setNodes(nodeList);

        e0->setId(nEdgeEl + i);
        int physicalTag = ng_mesh->GetFaceIndexOfSE(i);

        std::shared_ptr<MeshEntity> ent;
        if (face_id.find(physicalTag) != face_id.end())
            ent = copied_mesh->getEntityById(face_id[physicalTag]);
        else
            ent = copied_mesh->getEntityById(physicalTag);

        if (ent == nullptr)
        {
            std::string message = "Failed to copy mesh data. Please check the input geometry! (getEntityById returned nullptr)";
            LOG(ERROR) << message << physicalTag;

            throw cenos_exception(message);
        }

        e0->setEntity(ent);

        copied_mesh->addElement(e0);
    }

    //volume elements
    for (int i = 1; i <= nVolEl; i++)
    {
        int en[8]; //element nodes
        ElementTypes elType = ng_mesh->GetVolumeElement(i, en);

        std::shared_ptr<Element> e0 = Element::createElement(elType);

        std::vector<std::shared_ptr<Node>> nodeList;
        for (int j = 0; j < e0->getTypeNodeCount(); j++)
        {
            nodeList.push_back(copied_mesh->getNode(en[j]));
        }

        e0->setNodes(nodeList);

        e0->setId(i + nFaceEl + nEdgeEl);

        int physicalTag = ng_mesh->GetSolidIndexOfElement(i);

        std::shared_ptr<MeshEntity> ent;
        if (solid_id.find(physicalTag) != solid_id.end())
            ent = copied_mesh->getEntityById(solid_id[physicalTag]);
        else
            ent = copied_mesh->getEntityById(physicalTag); 
        

        if (ent == nullptr)
        {
            std::string message = "Failed to copy mesh data. Please check the input geometry! (getEntityById returned nullptr)";
            throw cenos_exception(message);
        }

        e0->setEntity(ent);

        copied_mesh->addElement(e0);
    }

    return copied_mesh;
}

// Make sure that circles are not meshed with just one segmnet, making them into lines
MeshParameters_ checkCurveSegments(TopoEntity_ ent, MeshParameters_ mp) {

    if (ent->getType() == ShapeType::EDGE) {
        Standard_Real first = 0, last = 0;
        Handle_Geom_Curve theCurve = BRep_Tool::Curve(TopoDS::Edge(ent->getShape()), first, last);
        if (theCurve->DynamicType() != STANDARD_TYPE(Geom_Line)) {

            double curveSpanAngle = last - first;
            // if curve spans more than 120 degrees (in radians), make the segments not span more than that
            if (curveSpanAngle > 2.0944) {
                GProp_GProps curveProp;
                BRepGProp::LinearProperties(ent->getShape(), curveProp);
                Standard_Real curveLength = curveProp.Mass();

                double oldSize = mp->getMaxH();
                double newSize = curveLength / (curveSpanAngle / 2.0944);
                if (newSize < oldSize) {
                    // We have to create a new pointer, otherwise it will affect the
                    // sizing of all other edges in that domain which have the same mp.
                    MeshParameters_ newMp = std::make_shared<MeshParameters>(*mp);
                    newMp->setMaxH(newSize);
                    newMp->setMinH(0.0);

                    std::cout << "Mesh size had to be reduced for " << ent->getName() << "(" << theCurve->DynamicType()->Name() << ")" << " from " << oldSize << " to " << newSize << std::endl;
                    LOG(INFO) << "Mesh size had to be reduced for " << ent->getName() << "(" << theCurve->DynamicType()->Name() << ")" << " from " << oldSize << " to " << newSize << std::endl;
                    return newMp;
                }
                
            }
        }
    }
    // All curves have good sizing, return the original mp
    return mp;
}

NetgenMesh_ MeshGenerator::generate1DMesh(std::shared_ptr<TopoEntity> ent)
{
    NetgenGeometryWrapper ng_geom(ent->getShape());

    NetgenMesh_ main_mesh = std::make_shared<NetgenMesh>();

    MeshParameters_ mp1 = checkCurveSegments(ent, meshing_data.getParameters(ent));
    auto mp = mp1;

    main_mesh->SetMeshSizing(&ng_geom, mp);
    main_mesh->CopyLocalH(main_ng_mesh);

    // if already meshed, return existing mesh
    if (meshing_data.isMeshed(ent))
    {
        main_mesh->MergeMesh(meshing_data.getMesh(ent));
        return main_mesh;
    }

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    NetgenLocalH* localh = meshing_data.getLocalH(ent);
    if (localh != nullptr)
    {
        main_mesh->RestrictSizeByLocalH(localh);
    }

    try
    {
        main_mesh->DivideEdges(ng_geom, *mp);
    }
    catch (std::exception e)
    {
        main_mesh->DeleteMesh();
        std::string message = "Could not create edge mesh. Please check the input geometry! (generate1DMesh(" + ent->getName() + ")";
        throw cenos_exception(message);
    }

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    meshing_data.setMesh(ent, main_mesh);
    return main_mesh;
}

NetgenMesh_ MeshGenerator::generate2DMesh(std::shared_ptr<TopoEntity> ent)
{
    TopoDS_Shape sh = ent->getShape();

    if (sh.IsNull())
    {
        std::string message = "Shape of entity " + ent->getName() + " is NULL!!!";
        throw cenos_exception(message);
    }
    NetgenGeometryWrapper ng_geom(sh);

    NetgenMesh_ main_mesh = std::make_shared< NetgenMesh>();

    auto mp = meshing_data.getParameters(ent);
    main_mesh->SetMeshSizing(&ng_geom, mp);
    main_mesh->CopyLocalH(main_ng_mesh);

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    // if already meshed, return existing mesh
    if (meshing_data.isMeshed(ent))
    {
        main_mesh->MergeMesh(meshing_data.getMesh(ent));
        return main_mesh;
    }

    TopTools_IndexedMapOfShape wmap;
    TopExp::MapShapes(sh, TopAbs_WIRE, wmap);
    std::map<int, int> edges_in_face; // map entity id to local edge id 
    int current_edge_index = 0;
    std::vector<int> added_entity_ids; //used to check if edge already added (e.g. seam edge)
    for (int w = 1; w <= wmap.Extent(); w++)
    {
        //loop all edges, mesh them and merge into current mesh (occ_mesh)
        TopTools_IndexedMapOfShape emap;
        TopExp::MapShapes(wmap.FindKey(w), TopAbs_EDGE, emap);

        for (TopExp_Explorer edge_explorer(wmap.FindKey(w), TopAbs_EDGE); edge_explorer.More(); edge_explorer.Next())
        {
            TopoDS_Shape current_edge = edge_explorer.Current();
            TopoDS_Edge edge = TopoDS::Edge(current_edge);

            if (BRep_Tool::Degenerated(edge))
                continue;

            
            TopoEntity_ edge_entity = meshing_data.getTopoEntity(current_edge);
            if (edge_entity == nullptr)
            {
                writeDebugFiles(wDir, current_edge, "failed_getTopoEntity_edge.brep");
                writeDebugFiles(wDir, ent->getShape(), "failed_mesh_ent.brep");
                main_mesh->DeleteMesh();
                throw cenos_exception("Could not getTopoEntity " + ent->getName() + " for 2D mesh. Please check input geometry!");
            }

            if (std::find(added_entity_ids.begin(), added_entity_ids.end(), edge_entity->getId()) != added_entity_ids.end())
                continue;

            double s0, s1;
            Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);

            NetgenMesh_ edge_mesh = generate1DMesh(edge_entity);

            bool is_reversed = false;
            if (current_edge.Orientation() != edge_entity->getShape().Orientation())
            {
                edge_mesh->ReverseSegments();
                is_reversed = true;
            }
            //this is due to recent changes in netgen. Now edges with reversed orientation should also have reversed segment orientation
            if (current_edge.Orientation() == TopAbs_Orientation::TopAbs_REVERSED)
            {
                edge_mesh->ReverseSegments();
                if (is_reversed)
                    is_reversed = false;
                else
                    is_reversed = true;
            }       

            if (stopper != nullptr and *stopper)
                throw cenos_exception("Stopped by user.");

            main_mesh->MergeMesh(edge_mesh, current_edge_index);
            current_edge_index++;
            added_entity_ids.push_back(edge_entity->getId());
            main_mesh->RestrictSizeByMesh(edge_mesh);

            //reverse back so that original mesh stored in mesh_data is not changed
            if (is_reversed)
                edge_mesh->ReverseSegments();

        }
    }

    // This function calls CalcSurfacesOfNode from netgen core.
    // Not sure what this one does, but it just does not work without it!
    // It has to be done if edges are copied into mesh instead of using GenerateEdgeMesh.
    // V.G. 19/Oct/2020
    main_mesh->CalculateSurfaceOfNode();
    //main_mesh->PrepareSurfaceMeshing();

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    NetgenLocalH* localh = meshing_data.getLocalH(ent);
    if (localh != nullptr)
    {
        main_mesh->RestrictSizeByLocalH(localh);
    }

    try
    {
        main_mesh->GenerateSurfaceMesh(ng_geom, *mp);
        
        //generate 2D boundary layer only if mesh is 2D
        if (mesh_mode == M_2D)
        {
            auto thicknesses = meshing_data.getParameters(ent)->getLayers();
            if (!thicknesses.empty())
            {
                double* th = &thicknesses[0];
                main_mesh->GenerateBoundaryLayer2D(1, th, thicknesses.size());
            }
        }
    }
    catch (std::exception e)
    {
        misc::createFolder(wDir + "/debug");
        std::string fn = wDir + "/debug/face_mesh.vol";
        main_mesh->SaveMesh(fn.c_str());
        ShapeAlgorithms::writeBrepFile(ent->getShape(), wDir + "/debug/failed_mesh_face.brep");
        main_mesh->DeleteMesh();
        std::string message = "Could not create face mesh. Please check the input geometry! (generate2DMesh(" + ent->getName() + ")";
        throw cenos_exception(message);
    }

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    meshing_data.setMesh(ent, main_mesh);
    return main_mesh;
}

NetgenMesh_ MeshGenerator::generate3DMesh(std::shared_ptr<TopoEntity> ent)
{
    TopoDS_Shape sh = ent->getShape();

    NetgenGeometryWrapper ng_geom(sh);

    NetgenMesh_ main_mesh = std::make_shared< NetgenMesh>();

    if (meshing_data.isMeshed(ent))
    {
        main_mesh->MergeMesh(meshing_data.getMesh(ent));
        return main_mesh;
    }

    main_mesh->SetGeometry(ng_geom);
    auto mp = meshing_data.getParameters(ent);
    main_mesh->SetMeshSizing(&ng_geom, mp);
    main_mesh->CopyLocalH(main_ng_mesh);

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    //loop over shells
    TopTools_IndexedMapOfShape shellmap;
    TopExp::MapShapes(sh, TopAbs_SHELL, shellmap);
    for (int shell_i = 1; shell_i <= shellmap.Extent(); shell_i++)
    {
        //loop all faces, mesh them and merge into current mesh (occ_mesh)
        TopTools_IndexedMapOfShape fmap;
        TopExp::MapShapes(shellmap.FindKey(shell_i), TopAbs_FACE, fmap);
        for (int f = 1; f <= fmap.Extent(); f++)
        {
            TopoDS_Shape current_face = fmap.FindKey(f);

            auto face_entity = meshing_data.getTopoEntity(current_face);
            if (face_entity == nullptr)
            {
                writeDebugFiles(wDir, current_face, "failed_getTopoEntity_face.brep");
                writeDebugFiles(wDir, ent->getShape(), "failed_mesh_ent.brep");
                main_mesh->DeleteMesh();
                throw cenos_exception("Could not getTopoEntity " + ent->getName() + " for 3D mesh. Please check input geometry!");
            }
            NetgenMesh_ face_mesh = generate2DMesh(face_entity);
                                 
            bool is_reversed = false;
            if (current_face.Orientation() != face_entity->getShape().Orientation() and current_face.Orientation() != TopAbs_Orientation::TopAbs_INTERNAL)
            {
                face_mesh->ReverseFaces();
                is_reversed = true;
            }

            if (stopper != nullptr and *stopper)
                throw cenos_exception("Stopped by user.");

            main_mesh->MergeMesh(face_mesh);
            main_mesh->CalculateSurfaceOfNode();

            if (is_reversed)
                face_mesh->ReverseFaces();
        }
    }

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    NetgenLocalH* localh = meshing_data.getLocalH(ent);
    if (localh != nullptr)
    {
        main_mesh->RestrictSizeByLocalH(localh);
    }

    try
    {
        main_mesh->SaveMesh("D:/test/solid_before.vol");

        main_mesh->GenerateVolumeMesh(*mp);

        auto thicknesses = meshing_data.getParameters(ent)->getLayers();
        if (!thicknesses.empty())
        {
            //get surface id array
            TopTools_IndexedMapOfShape fmap;
            TopExp::MapShapes(sh, TopAbs_FACE, fmap);
            int* f_ids = new int[fmap.Extent()];
            for (int f = 1; f <= fmap.Extent(); f++)
            {
                f_ids[f - 1] = f;
            }

            double* th = &thicknesses[0];

            main_mesh->GenerateBoundaryLayers(f_ids, fmap.Extent(), th, thicknesses.size());
            delete[] f_ids;
        }
    }

    catch (std::exception e)
    {
        writeDebugFiles(wDir, main_mesh, "failed_occ_mesh.vol");
        writeDebugFiles(wDir, ent->getShape(), "failed_occ_mesh_entity.vol");
        main_mesh->DeleteMesh();
        std::string message = "Could not create volume mesh. Please check the input geometry! (generate3DMesh(" + ent->getName() + ")";
        throw cenos_exception(message);
    }

    if (stopper != nullptr and *stopper)
        throw cenos_exception("Stopped by user.");

    meshing_data.setMesh(ent, main_mesh);
    return main_mesh;

}


void MeshGenerator::validateInput()
{
    int edges = 0, faces = 0, solids = 0;
    for (auto ent : entities)
    {
        if (ent->getType() == ShapeType::EDGE)
            edges++;
        else if (ent->getType() == ShapeType::FACE)
            faces++;
        else if (ent->getType() == ShapeType::SOLID)
            solids++;
    }

    TopTools_IndexedMapOfShape emap;
    TopExp::MapShapes(partition, TopAbs_EDGE, emap);
    int nr_edges = 0;
    for (int e = 1; e <= emap.Extent(); e++)
    {
        TopoDS_Edge edge = TopoDS::Edge(emap.FindKey(e));
        if (!BRep_Tool::Degenerated(edge))
            nr_edges++;
    }

    if (nr_edges != edges)
    {
        std::string message = "Meshing failed! Check input geometry! Number of edges in passed geometry " + std::to_string(nr_edges) + " is not equal to number of TopoEntities of EDGE type (" + std::to_string(edges) + ").";
        throw cenos_exception(message);
    }


    TopTools_IndexedMapOfShape fmap;
    TopExp::MapShapes(partition, TopAbs_FACE, fmap);
    if (fmap.Extent() != faces)
    {
        std::string message = "Meshing failed! Check input geometry! Number of faces in passed geometry " + std::to_string(fmap.Extent()) + " is not equal to number of TopoEntities of FACE type (" + std::to_string(faces) + ").";
        throw cenos_exception(message);
    }


    TopTools_IndexedMapOfShape smap;
    TopExp::MapShapes(partition, TopAbs_SOLID, smap);
    if (smap.Extent() != solids)
    {
        std::string message = "Meshing failed! Check input geometry! Number of solids in passed geometry " + std::to_string(smap.Extent()) + " is not equal to number of TopoEntities of SOLID type (" + std::to_string(solids) + ").";
        throw cenos_exception(message);
    }


    //check iternal faces
    for (auto s : entities)
    {
        if (s->getType() == ShapeType::SOLID)
        {
            if(ShapeAlgorithms::hasInternalFaces(s->getShape()))
            {
                std::string message = "Currently MeshGenerator does not support internal faces ";
                throw cenos_exception(message);
            }
        }
    }
}




void MeshGenerator::copyTransformedMeshData()
{
    if (entity_transform.size() == 0)
        return;

    // ****************************************************
// Copying mesh data for moved and non-modified parts
// ****************************************************

    for (auto trsf : entity_transform)
    {
        auto ent = trsf.first;

        bool isDegenerated = false;
        if (ent->getShape().ShapeType() == TopAbs_ShapeEnum::TopAbs_EDGE) {
            if (BRep_Tool::Degenerated(TopoDS::Edge(ent->getShape()))) {
                isDegenerated = true;
            }
        }

        if (not isDegenerated) {
            // CNg_OCC_ShapeToGeometry will crash in case of degenerated edges. CP-1490
            NetgenMesh_ entity_mesh = mesh->getNetgenMesh(ent->getName(), trsf.second);
            NetgenGeometryWrapper entity_geom = ent->getShape();
            entity_mesh->SetGeometry(entity_geom);
            meshing_data.setMesh(ent, entity_mesh);
        }
    }

    // ****************************************************
    // Assigning mesh sizing for modified parts
    // ****************************************************

    for (auto ent : entities)
    {
        if (!ent->isEnabled())
            continue;

        if (entity_transform.find(ent) != entity_transform.end())
            continue;

        NetgenMesh_ entity_mesh = mesh->getNetgenMesh(ent->getName());

        MeshParameters_ mp = std::make_shared<MeshParameters>();
        double maxh = entity_mesh->GetMaxH();
        mp->setMaxH(maxh);
        mp->setMinH(0);
        mp->setGrading(0.4);
        meshing_data.setParameters(ent, mp);
        meshing_data.setLocalH(ent, entity_mesh->GetLocalH().get());
    }
}