#define NOMINMAX
#include <mystdlib.h>
#include <WinSock2.h>
#include <Windows.h>
#include <meshing.hpp>
#include <occgeom.hpp>
#include <Standard_Version.hxx>

#include <nginterface.h>
#include <myadt.hpp>

#include "NetgenWrapper.h"
#include "misc/miscFunctions.hpp"
#include "cenos_exception.h"
#include <TopoDS_Edge.hxx>
#include <TopExp.hxx>

std::ofstream logstream;

namespace netgen
{
    DLL_HEADER extern OCCParameters occparam;

    DLL_HEADER extern MeshingParameters mparam;
}

using namespace netgen;

void DivideEdge(std::shared_ptr<NgMesh> mesh, TopoDS_Edge& edge, netgen::NgArray<netgen::MeshPoint>& ps, netgen::NgArray<double>& params);

void InitializeNetgen()
{
    std::string filename = misc::getLogPrefix() + ".cenos-meshing.log";
    logstream.open(filename, std::ofstream::out | std::ofstream::app);
    if (logstream.is_open())
    {
        netgen::mycout = &logstream;
    }
    else
    {
        netgen::mycout = &cout;
        netgen::myerr = &cerr;
        std::cerr << "Could not open meshing.log file. Meshing log will be redrected to cout" << std::endl;
    }
}

void FinalizeNetgen()
{
    logstream.close();
}



void TransferSizing(MeshParameters mp)
{
    mparam.uselocalh = 1;
   // mparam.closeedgefac = 2.0;			    // Factor to use for refinement at close edges

    mparam.maxh = mp.getMaxH() * mp.getDensityFactor();
    mparam.minh = mp.getMinH() * mp.getDensityFactor();

    mparam.grading = mp.getGrading();
    mparam.curvaturesafety = 2.0;
    mparam.segmentsperedge = 2.0;

    mparam.secondorder = 0;
    mparam.quad = 0;

    mparam.meshsizefilename = "";
    mparam.optsteps2d = 3;
    mparam.optsteps3d = 3;

    mparam.inverttets = 0;
    mparam.inverttrigs = 0;

    mparam.checkoverlap = 1;
    mparam.checkoverlappingboundary = 1;
}


NetgenMesh::NetgenMesh()
{
	mesh = std::make_shared<NgMesh>();
	//mesh->AddFaceDescriptor(netgen::FaceDescriptor(1, 1, 0, 1));
}

void NetgenMesh::ReverseSegments()
{
	for (netgen::SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
	{
		swap((*mesh)[si][0], (*mesh)[si][1]);
		swap((*mesh)[si].epgeominfo[0].dist, (*mesh)[si].epgeominfo[1].dist);
		swap((*mesh)[si].epgeominfo[0].edgenr, (*mesh)[si].epgeominfo[1].edgenr);
		swap((*mesh)[si].epgeominfo[0].u, (*mesh)[si].epgeominfo[1].u);
		swap((*mesh)[si].epgeominfo[0].v, (*mesh)[si].epgeominfo[1].v);
	}
}

void NetgenMesh::MergeMesh(NetgenMesh_ orig_mesh, int index)
{
    //Geometries
    const netgen::OCCGeometry& origin_geom = *dynamic_pointer_cast<netgen::OCCGeometry>(orig_mesh->mesh->GetGeometry());
    const netgen::OCCGeometry& this_geom = *dynamic_pointer_cast<netgen::OCCGeometry> (this->mesh->GetGeometry());
    double eps = 1e-6 * this_geom.GetBoundingBox().Diam();

    //if origin is edge
    int edge_index = 0;
    int solid_index = 0;
    int origin_mesh_dimension = 0;

    if (origin_geom.emap.Extent() == 1 && origin_geom.fmap.Extent() == 0 && origin_geom.somap.Extent() == 0)
    {
        if (index == -1)
            edge_index = this_geom.emap.FindIndex(origin_geom.emap.FindKey(1)) - 1;
        else
            edge_index = index;
        origin_mesh_dimension = 1;
    }
    else if (origin_geom.fmap.Extent() == 1 && origin_geom.somap.Extent() == 0) // origin is face
    {
        if (index == -1)
            this->mesh->AddFaceDescriptor(netgen::FaceDescriptor(this->mesh->GetNFD() + 1, 1, 0, 1));
        else
            this->mesh->AddFaceDescriptor(netgen::FaceDescriptor(index, 1, 0, 1));
        origin_mesh_dimension = 2;
    }
    else if (origin_geom.somap.Extent() == 1)
    {
        if (index == -1)
            solid_index = this_geom.somap.FindIndex(origin_geom.somap.FindKey(1));
        else
            solid_index = index;
        origin_mesh_dimension = 3;
    }
    else // not supported
    {
        std::string message = "Ng_MergeMesh() supports only single edge or single face or single solid as origin geometry. Passed Edges: " + std::to_string(origin_geom.emap.Extent())
            + ", Faces: " + std::to_string(origin_geom.fmap.Extent()) + ", Solids: " + std::to_string(origin_geom.somap.Extent()) + ".";
        std::cout << message << std::endl;
        throw(message);
    }

    // map to hold node ids.
    // Example:
    // node 1 will correspond to node 10, if there are already 9 nodes in destination mesh
    // node 3 will correspond to node 7, if they are equal (say, edge node)
    std::map<netgen::PointIndex, netgen::PointIndex> node_ids;

    int nNodes = orig_mesh->mesh->GetNP();
    for (netgen::PointIndex opi : orig_mesh->mesh->Points().Range())
    {
        const netgen::Point3d& p = (*orig_mesh->mesh)[opi];
        bool exists = false;
        for (netgen::PointIndex dpi : this->mesh->Points().Range())
        {
            if (netgen::Dist2((*this->mesh)[dpi], p) < eps * eps)
            {
                exists = true;
                node_ids[opi] = dpi;
                break;
            }
        }

        if (!exists)
        {
            this->mesh->AddPoint(p);
            node_ids[opi] = netgen::PointIndex(this->mesh->GetNP());
        }
    }

    int nEdgeEl = orig_mesh->mesh->GetNSeg();
    int nFaceEl = orig_mesh->mesh->GetNSE();
    int nVolEl = orig_mesh->mesh->GetNE();

    //edge elements
    if (origin_mesh_dimension == 1)
    {
        for (int i = 1; i <= nEdgeEl; i++)
        {
            netgen::Segment seg = orig_mesh->mesh->LineSegment(i);

            seg[0] = node_ids[seg[0]];
            seg[1] = node_ids[seg[1]];
            seg.epgeominfo[0].edgenr = edge_index;
            seg.epgeominfo[1].edgenr = edge_index;
            seg.edgenr = edge_index + 1;
            seg.si = 1;
            seg.cd2i = -1;

            this->mesh->AddSegment(seg);
        }
    }

    //surface elements
    if (origin_mesh_dimension == 2)
    {
        for (int i = 1; i <= nFaceEl; i++)
        {
            //TODO check element type!!!
            netgen::Element2d el = orig_mesh->mesh->SurfaceElement(i);
            el.PNum(1) = node_ids[el.PNum(1)];
            el.PNum(2) = node_ids[el.PNum(2)];
            el.PNum(3) = node_ids[el.PNum(3)];
            if (el.GetType() == netgen::QUAD)
                el.PNum(4) = node_ids[el.PNum(4)];

            el.SetIndex(this->mesh->GetNFD());
            this->mesh->AddSurfaceElement(el);
        }
    } 


    //volume elements
    if (origin_mesh_dimension == 3)
    {
        for (int i = 1; i <= nVolEl; i++)
        {
            netgen::Element el = orig_mesh->mesh->VolumeElement(i);
            for (int n = 1; n <= el.GetNP(); n++)
                el.PNum(n) = node_ids[el.PNum(n)];
            el.SetIndex(solid_index);
            this->mesh->AddVolumeElement(el);
        }
    }
}

void NetgenMesh::ReverseFaces()
{
    for (netgen::SurfaceElementIndex si = 0; si < mesh->GetNSE(); si++)
    {
        (*mesh)[si].Invert();
    }
}

void NetgenMesh::DeleteMesh()
{
    if (mesh != NULL)
    {
        // Delete the Mesh structures
        mesh->DeleteMesh();

        // Set mesh new ptr
        mesh = std::make_shared<NgMesh>();
    }
}

void NetgenMesh::SetElementRefinement(int el_index, bool flag)
{
    if (mesh->GetDimension() == 3)
    {
        mesh->VolumeElement(el_index).SetRefinementFlag(flag);
    }
    else
    {
        mesh->SurfaceElement(el_index).SetRefinementFlag(flag);
    }
}

void NetgenMesh::Refine()
{
    netgen::BisectionOptions biopt;
    biopt.usemarkedelements = 1;
    biopt.refine_p = 0;
    biopt.refine_hp = 0;

    // Set refinement false for surface elements.
    // TO DO what to do for 2D meshes?
    for (int i = 1; i <= mesh->GetNSE(); i++)
        mesh->SurfaceElement(i).SetRefinementFlag(false);

    mesh->GetGeometry()->GetRefinement().Bisect(*mesh, biopt);
    mesh->UpdateTopology();
    mesh->GetCurvedElements().SetIsHighOrder(false);
}


void NetgenMesh::GetPoint(int index, double* pt)
{
    const netgen::Point3d& p = mesh->Point(index);
    pt[0] = p.X();
    pt[1] = p.Y();
    pt[2] = p.Z();
}

void NetgenMesh::GetSegment(int index, int* point_indices)
{
    const netgen::Segment& seg = mesh->LineSegment(index);
    point_indices[0] = seg[0];
    point_indices[1] = seg[1];
}

void NetgenMesh::SetGeometry(NetgenGeometryWrapper ng_geom)
{
    mesh->SetGeometry(shared_ptr<netgen::OCCGeometry>(ng_geom.occ_shape));
}

double NetgenMesh::GetMaxH()
{
    mesh->CalcLocalHFromPointDistances(0.4);

    double max_h = 0;
    for (netgen::PointIndex opi : mesh->Points().Range())
    {
        const netgen::Point3d& p = (*mesh)[opi];
        double point_h = mesh->GetH(p);
        if (point_h > max_h)
            max_h = point_h;
    }
    return max_h;
}

NetgenMesh::~NetgenMesh()
{
   // DeleteMesh();
}

void NetgenMesh::AddPoint(double* pt)
{
    mesh->AddPoint(netgen::Point3d(pt[0], pt[1], pt[2]));
}

void NetgenMesh::AddSegment(int pi1, int pi2, int edgeIndex)
{
    const netgen::Point3d& p1 = mesh->Point(pi1);
    const netgen::Point3d& p2 = mesh->Point(pi2);

    netgen::Segment seg;
    seg[0] = pi1;
    seg[1] = pi2;

    seg.epgeominfo[0].edgenr = edgeIndex;
    seg.epgeominfo[1].edgenr = edgeIndex;

    seg.epgeominfo[0].u = p1.X();
    seg.epgeominfo[0].v = p1.Y();
    seg.epgeominfo[1].u = p2.X();
    seg.epgeominfo[1].v = p2.X();

 //   seg.si = edgeIndex + 1;
    seg.edgenr = edgeIndex + 1;

    mesh->AddSegment(seg);
}

void NetgenMesh::AddSurfaceElement(ElementTypes et, int* pi, int surfIndx)
{
    Element2d el(3);
    if (et == ElementTypes::TRI)
        el = Element2d(3);
    else if (et == ElementTypes::QUAD)
        el = Element2d(4);

    el.SetIndex(surfIndx);
    el.PNum(1) = pi[0];
    el.PNum(2) = pi[1];
    el.PNum(3) = pi[2];

    //add fourth node if it is quad
    if (et == ElementTypes::QUAD)
        el.PNum(4) = pi[3];

    mesh->AddSurfaceElement(el);
}

void NetgenMesh::AddVolumeElement(ElementTypes et, int* pi, int volIndx)
{
    int nodeCount = 4;
    if (et == ElementTypes::PRISM)
    {
        nodeCount = 6;
    }

    netgen::Element el(nodeCount);
    el.SetIndex(volIndx);
    el.PNum(1) = pi[0];
    el.PNum(2) = pi[1];
    el.PNum(3) = pi[2];
    el.PNum(4) = pi[3];
    if (et == ElementTypes::PRISM)
    {
        el.PNum(5) = pi[4];
        el.PNum(6) = pi[5];
    }

    mesh->AddVolumeElement(el);
}



void NetgenMesh::AddFaceDescriptor(int faceId, int dominId, int domOutId, int TLOface)
{
    mesh->AddFaceDescriptor(netgen::FaceDescriptor(faceId, dominId, domOutId, TLOface));
}



void NetgenMesh::SetMeshSizing(NetgenGeometryWrapper* ng_geom, MeshParameters_ mp)
{
    OCCGeometry* occgeom = (ng_geom->occ_shape).get();

    mesh->SetGeometry(ng_geom->occ_shape);

    mesh->geomtype = NgMesh::GEOM_OCC;

    TransferSizing(*mp);

    //clean up mesh
    mesh->DeleteMesh();

    OCCSetLocalMeshSize(*(ng_geom->occ_shape), *mesh, mparam, occparam);
}


std::shared_ptr<NetgenLocalH> NetgenMesh::GetLocalH()
{
    mesh->CalcLocalHFromPointDistances(0.4);
    return mesh->GetLocalH();
}


void NetgenMesh::PrepareSurfaceMeshing()
{
    const netgen::OCCGeometry& occgeom = *dynamic_pointer_cast<netgen::OCCGeometry>(mesh->GetGeometry());

    TopoDS_Face face = TopoDS::Face(occgeom.fmap.FindKey(1));
    int uvc = 0;
    for (int e = 1; e <= occgeom.emap.Extent(); e++)
    {
        TopoDS_Edge edge = TopoDS::Edge(occgeom.emap.FindKey(e));
        Handle(Geom2d_Curve) cof;
        double s0, s1;
        cof = BRep_Tool::CurveOnSurface(edge, face, s0, s1);

        for (netgen::SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
        {
            if ((*mesh)[si].epgeominfo[1].edgenr == e)
            {
                gp_Pnt2d p2d;
                uvc++;
                p2d = cof->Value((*mesh)[si].epgeominfo[0].dist);

                (*mesh)[si].epgeominfo[0].u = p2d.X();

                (*mesh)[si].epgeominfo[0].v = p2d.Y();

                p2d = cof->Value((*mesh)[si].epgeominfo[1].dist);
                (*mesh)[si].epgeominfo[1].u = p2d.X();
                (*mesh)[si].epgeominfo[1].v = p2d.Y();
            }

        }

    }
}

void NetgenMesh::CalculateSurfaceOfNode()
{
    mesh->CalcSurfacesOfNode();
}

void NetgenMesh::RestrictSizeByLocalH(NetgenLocalH* localh)
{
    for (netgen::PointIndex opi : mesh->Points().Range())
    {
        const netgen::Point3d& p = (*mesh)[opi];
        mesh->RestrictLocalH(p, localh->GetH(p));
    }
}

void NetgenMesh::RestrictSizeByMesh(NetgenMesh_ ng_mesh)
{
    std::shared_ptr<NgMesh> origin_mesh = ng_mesh->mesh;

    int nNodes = origin_mesh->GetNP();
    for (netgen::PointIndex opi : origin_mesh->Points().Range())
    {
        const netgen::Point3d& p = (*origin_mesh)[opi];
        mesh->RestrictLocalH(p, origin_mesh->GetH(p));
    }
}

void NetgenMesh::CopyLocalH(NetgenMesh_ ng_mesh)
{
    std::shared_ptr<NgMesh> origin_mesh = ng_mesh->mesh;
    mesh->SetLocalH(origin_mesh->GetLocalH());
}

void NetgenMesh::SaveMesh(const char* filename)
{
    mesh->Save(filename);
}

int NetgenMesh::GetNP()
{
    return mesh->GetNP();
}

int NetgenMesh::GetNSeg()
{
    return mesh->GetNSeg();
}

int NetgenMesh::GetNSE()
{
    return mesh->GetNSE();
}

int NetgenMesh::GetNE()
{
    return mesh->GetNE();
}

int NetgenMesh::GetEdgeIndexOfSegment(int i)
{
    return mesh->LineSegment(i).epgeominfo[0].edgenr;
}

int NetgenMesh::GetFaceIndexOfSE(int i)
{
    int fd_id =mesh->SurfaceElement(i).GetIndex();
    return mesh->GetFaceDescriptor(fd_id).SurfNr();
}


int NetgenMesh::GetSolidIndexOfElement(int i)
{
    return mesh->VolumeElement(i).GetIndex();
}

ElementTypes NetgenMesh::GetSurfaceElement(int i, int* pi)
{
    const netgen::Element2d& el = mesh->SurfaceElement(i);
    for (int i = 1; i <= el.GetNP(); i++)
        pi[i - 1] = el.PNum(i);
    ElementTypes et;
    switch (el.GetNP())
    {
    case 3: et = ElementTypes::TRI; break;
    case 4: et = ElementTypes::QUAD; break;
    case 6:
        switch (el.GetNV())
        {
        case 3: et = ElementTypes::TRI_2; break;
        default:
            et = ElementTypes::TRI_2; break;
        }
        break;
    case 8: et = ElementTypes::QUAD_2_8N; break;
    default:
        et = ElementTypes::TRI; break;
    }
    return et;
}

ElementTypes NetgenMesh::GetVolumeElement(int i, int* pi)
{
    const netgen::Element& el = mesh->VolumeElement(i);
    for (int i = 1; i <= el.GetNP(); i++)
        pi[i - 1] = el.PNum(i);
    ElementTypes et;
    switch (el.GetNP())
    {
    case 4: et = ElementTypes::TETRA; break;
    case 5: et = ElementTypes::PYRA; break;
    case 6: et = ElementTypes::PRISM; break;
    case 10: et = ElementTypes::TETRA_2; break;
    default:
        et = ElementTypes::TETRA; break;
    }
    return et;
}

void NetgenMesh::GenerateEdgeMesh(NetgenGeometryWrapper ng_geom, MeshParameters mp)
{
    mesh->SetGeometry(ng_geom.occ_shape);

    TransferSizing(mp);
    ng_geom.occ_shape->FindEdges(*mesh, mparam);


}

void NetgenMesh::DivideEdges(NetgenGeometryWrapper ng_geom, MeshParameters mp)
{
    mesh->SetGeometry(ng_geom.occ_shape);
    TransferSizing(mp);
    double eps = 1e-6 * ng_geom.occ_shape->GetBoundingBox().Diam();
    
    // add vertices to mesh as points
    for (int i = 1; i <= ng_geom.occ_shape->vmap.Extent(); i++)
    {
        gp_Pnt pnt = BRep_Tool::Pnt(TopoDS::Vertex(ng_geom.occ_shape->vmap(i)));
        netgen::MeshPoint mpnt(netgen::Point<3>(pnt.X(), pnt.Y(), pnt.Z()));
        
        bool exists = 0;
        for (netgen::PointIndex pi : mesh->Points().Range())
        {
            if (netgen::Dist2((*mesh)[pi], netgen::Point<3>(mpnt)) < eps * eps)
            {
                exists = true;
                break;
            }
        }
        if (!exists)
            mesh->AddPoint(mpnt);
    }


    netgen::PointIndex first_ep = mesh->Points().Range().end();
    auto vertexrange = mesh->Points().Range();
    for (int edge_i = 1; edge_i <= ng_geom.occ_shape->emap.Extent(); edge_i++)
    {
        TopoDS_Edge edge = TopoDS::Edge(ng_geom.occ_shape->emap(edge_i));
        if (BRep_Tool::Degenerated(edge))
        {
            //(*testout) << "ignoring degenerated edge" << endl;
            continue;
        }

        int geomedgenr = ng_geom.occ_shape->emap.FindIndex(edge);
        if (geomedgenr < 1) continue;

        if (ng_geom.occ_shape->vmap.FindIndex(TopExp::FirstVertex(edge)) ==
            ng_geom.occ_shape->vmap.FindIndex(TopExp::LastVertex(edge)))
        {
            GProp_GProps system;
            BRepGProp::LinearProperties(edge, system);

            if (system.Mass() < eps)
            {
                cout << "ignoring edge " << ng_geom.occ_shape->emap.FindIndex(edge)
                    << ". closed edge with length < " << eps << endl;
                continue;
            }
        }

        netgen::NgArray <netgen::MeshPoint> mp;
        netgen::NgArray <double> params;

        DivideEdge(mesh, edge, mp, params);
        netgen::NgArray<netgen::PointIndex> pnums(mp.Size() + 2);
        gp_Pnt p1 = BRep_Tool::Pnt(TopExp::FirstVertex(edge));
        netgen::Point<3> fp = netgen::Point<3>(p1.X(), p1.Y(), p1.Z());

        gp_Pnt p2 = BRep_Tool::Pnt(TopExp::LastVertex(edge));
        netgen::Point<3> lp = netgen::Point<3>(p2.X(), p2.Y(), p2.Z());

        pnums[0] = netgen::PointIndex::INVALID;
        pnums.Last() = netgen::PointIndex::INVALID;
        for (netgen::PointIndex pi : vertexrange)
        {
            if (Dist2((*mesh)[pi], fp) < eps * eps) pnums[0] = pi;
            if (Dist2((*mesh)[pi], lp) < eps * eps) pnums.Last() = pi;
        }

        for (size_t i = 1; i <= mp.Size(); i++)
        {
            bool exists = 0;
            for (netgen::PointIndex pi : vertexrange)
                if (((*mesh)[pi] - netgen::Point<3>(mp[i - 1])).Length() < eps)
                {
                    exists = true;
                    pnums[i] = pi;
                    break;
                }

            if (!exists)
                pnums[i] = mesh->AddPoint(mp[i - 1]);
        }

        (*netgen::mycout) << "NP = " << mesh->GetNP() << endl;

        for (size_t i = 1; i <= mp.Size() + 1; i++)
        {
            netgen::Segment seg;

            seg[0] = pnums[i - 1];
            seg[1] = pnums[i];
            seg.edgenr = geomedgenr + 1;
            seg.si = 1;
            seg.cd2i = -1;
            seg.epgeominfo[0].dist = params[i - 1];
            seg.epgeominfo[1].dist = params[i];
            seg.epgeominfo[0].edgenr = geomedgenr;
            seg.epgeominfo[1].edgenr = geomedgenr;

            if (edge.Orientation() == TopAbs_REVERSED)
            {
                swap(seg[0], seg[1]);
                swap(seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
                swap(seg.epgeominfo[0].u, seg.epgeominfo[1].u);
                swap(seg.epgeominfo[0].v, seg.epgeominfo[1].v);
            }

            mesh->AddSegment(seg);
        }
    }


    if ((mesh->GetNP()) && (mesh->GetNSeg()))
    {
        return;
    }
    else
    {
        throw(cenos_exception("No points or segments in 1D mesh"));
    }
}


void DivideEdge(std::shared_ptr<NgMesh> mesh, TopoDS_Edge& edge, netgen::NgArray<netgen::MeshPoint>& ps, netgen::NgArray<double>& params)
{
    double s0, s1;
    int nsubedges = 1;
    gp_Pnt pnt, oldpnt;
    const int DIVIDEEDGESECTIONS = 10000;
    double svalue[DIVIDEEDGESECTIONS];

    GProp_GProps system;
    BRepGProp::LinearProperties(edge, system);
    double L = system.Mass();

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

    double hvalue[DIVIDEEDGESECTIONS + 1];
    hvalue[0] = 0;
    pnt = c->Value(s0);

    int tmpVal = (int)(DIVIDEEDGESECTIONS);

    for (auto i : Range(DIVIDEEDGESECTIONS))
    {
        oldpnt = pnt;
        pnt = c->Value(s0 + double(i+1) / double(DIVIDEEDGESECTIONS) * (s1 - s0));
        auto pt3 = netgen::Point3d(pnt.X(), pnt.Y(), pnt.Z());
        hvalue[i+1] = hvalue[i ] +  1.0 / mesh->GetH(pt3) * pnt.Distance(oldpnt);
    }

    //  nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS]));
    nsubedges = max(1, int(floor(hvalue[DIVIDEEDGESECTIONS] + 0.5)));

    ps.SetSize(nsubedges - 1);
    params.SetSize(nsubedges + 1);

    int i = 1;
    int i1 = 0;
    do
    {
        if (hvalue[i1] / hvalue[DIVIDEEDGESECTIONS] * nsubedges >= i)
        {
            params[i] = (double(i1) / DIVIDEEDGESECTIONS);
            pnt = c->Value(s0 + params[i] * (s1 - s0));
            ps[i - 1] = netgen::MeshPoint(netgen::Point3d(pnt.X(), pnt.Y(), pnt.Z()));
            i++;
        }
        i1++;
        if (i1 > DIVIDEEDGESECTIONS)
        {
            nsubedges = i;
            ps.SetSize(nsubedges - 1);
            params.SetSize(nsubedges + 1);
            cout << "divide edge: local h too small" << endl;
        }
    } while (i < nsubedges);

    params[0] = 0;
    params[nsubedges] = 1.0;

    if (params[nsubedges] <= params[nsubedges - 1])
    {
        cout << "CORRECTED" << endl;
        ps.SetSize(nsubedges - 2);
        params.SetSize(nsubedges);
        params[nsubedges] = 1.0;
    }
}

void NetgenMesh::GenerateSurfaceMesh(NetgenGeometryWrapper ng_geom, MeshParameters mp)
{
    int numpoints = 0;
    (*netgen::mycout) << "Starting GenerateSurfaceMesh" << endl << flush;

    // if not face descriptors, add one
    mesh->AddFaceDescriptor(FaceDescriptor(1, 1, 0, 1));

    if (!mesh->GetNFD())
        mesh->AddFaceDescriptor(netgen::FaceDescriptor(1, 1, 0, 1));

    mesh->SetGeometry(ng_geom.occ_shape);
    (*netgen::mycout) << "Geometry is set" << endl << flush;

    // Set the internal meshing parameters structure from the nglib meshing 
    // parameters structure
    TransferSizing(mp);
    numpoints = mesh->GetNP();

    (*netgen::mycout) << "Parameters are transferred " << numpoints << endl << flush;
    

    // Initially set up only for surface meshing without any optimisation
    int perfstepsend = netgen::MESHCONST_MESHSURFACE;

    //optimize surface (we have not set a way to disable this)
    perfstepsend = netgen::MESHCONST_OPTSURFACE;


    (*netgen::mycout) << "Start mesh surface" << endl << flush;
    try
    {
        ng_geom.occ_shape->MeshSurface(*mesh, mparam);
        ng_geom.occ_shape->OptimizeSurface(*mesh, mparam);

    }
    catch (netgen::NgException e)
    {
        std::string msg = "Netgen Exception: " + std::string(e.what());
        throw(cenos_exception(msg));
    }

    (*netgen::mycout) << "Done meshing surface" << endl << flush;;

    mesh->CalcSurfacesOfNode();

    (*netgen::mycout) << mesh->GetNP() << endl << flush;
    (*netgen::mycout) << numpoints << endl << flush;
    (*netgen::mycout) << mesh->GetNSE() << endl << flush;

 // if required to make second order, use
  //      me->GetGeometry()->GetRefinement().MakeSecondOrder(*me);

    if (mesh->GetNP() < numpoints)
    {
        std::string msg = "Netgen Exception: nr of points smaller than initial.";
        std::cout << "Throwing " << std::endl;
        throw(cenos_exception(msg));
    }

    if (mesh->GetNSE() <= 0)
    {
        std::string msg = "Netgen Exception: no surface elements generated.";
        std::cout << "Throwing " << std::endl;

        throw(cenos_exception(msg));
    }

}

void NetgenMesh::GenerateBoundaryLayer2D(int dom_nr, double* heights_arr, int heights_count)
{

    netgen::Array<double> thicknesses;
    for (int i = 0; i < heights_count; i++)
    {
        thicknesses.Append(heights_arr[i]);
    }

    try
    {
        netgen::GenerateBoundaryLayer2(*mesh, dom_nr, thicknesses, false);
    }
    catch (netgen::NgException e)
    {
        std::string msg = "Netgen Exception: " + std::string(e.what());
        throw(cenos_exception(msg));
    }
}

void NetgenMesh::GenerateBoundaryLayers(int* surfid_arr, int surfid_count, double* heights_arr, int heights_count)
{
    netgen::BoundaryLayerParameters blp;
    for (int i = 0; i < surfid_count; i++) { blp.surfid.Append(surfid_arr[i]); }
    for (int i = 0; i < heights_count; i++) { blp.heights.Append(heights_arr[i]); }
    blp.new_mat = "BoundaryLayer";
    blp.domains.SetSize(2);
    blp.domains.Clear();
    blp.domains.SetBit(1);
    blp.outside = false;
    blp.grow_edges = true;
    std::cout << "Preparing BL done" << std::endl;
    netgen::GenerateBoundaryLayer(*mesh, blp);
    std::cout << "Meshing BL done" << std::endl;
}


void NetgenMesh::GenerateVolumeMesh(MeshParameters mp)
{
    TransferSizing(mp);

    mesh->CalcLocalH(mparam.grading);

    MeshVolume(mparam, *mesh);
    RemoveIllegalElements(*mesh);
    OptimizeVolume(mparam, *mesh);
}


NetgenGeometryWrapper::NetgenGeometryWrapper(TopoDS_Shape* shape)
{
	occ_shape = std::make_shared < NetgenOCCShape>();
	occ_shape->shape = *shape;
	occ_shape->changed = 1;
	occ_shape->BuildFMap();
	occ_shape->CalcBoundingBox();
}

NetgenGeometryWrapper::NetgenGeometryWrapper(TopoDS_Shape shape)
{
	occ_shape = std::make_shared < NetgenOCCShape>();
	occ_shape->shape = shape;
	occ_shape->changed = 1;
	occ_shape->BuildFMap();
	occ_shape->CalcBoundingBox();
}

TopTools_IndexedMapOfShape NetgenGeometryWrapper::GetSolidMap()
{
    return occ_shape->somap;
}

TopTools_IndexedMapOfShape NetgenGeometryWrapper::GetFaceMap()
{
    return occ_shape->fmap;
}

TopTools_IndexedMapOfShape NetgenGeometryWrapper::GetEdgeMap()
{
    return occ_shape->emap;
}