#include "cnglib.h"
#include <interface/writeuser.hpp>

#include <map>
#include "occgeom.hpp"
#include <linalg.hpp>
#include <csg.hpp>
#include <stlgeom.hpp>
#include <geometry2d.hpp>
#include <meshing.hpp>
#include <BRep_Builder.hxx>
//#include <../visualization/visual.hpp>

namespace netgen
{
    DLL_HEADER extern OCCParameters occparam;

    DLL_HEADER extern MeshingParameters mparam;
}

using namespace netgen;

inline void NOOP_Deleter(void*) { ; }


// Mapping of entities from Netgen definitions to GMSH definitions
enum GMSH_ELEMENTS {
    GMSH_LINE = 1, GMSH_TRIG = 2, GMSH_TRIG6 = 9,
    GMSH_QUAD = 3, GMSH_PRISM = 6, GMSH_QUAD8 = 16,
    GMSH_TET = 4, GMSH_TET10 = 11
};
const int triGmsh[7] = { 0,1,2,3,6,4,5 };
const int quadGmsh[9] = { 0,1,2,3,4,5,8,6,7 };
const int tetGmsh[11] = { 0,1,2,3,4,5,8,6,7,10,9 };
const int prismGmsh[6] = { 0,1,2,3,4,5 };
#define DIVIDEEDGESECTIONS 10000   // better solution to come soon

OCCGeometry* CreateOCCGeometryFromTopoDS(Cenos_TopoDS_Shape* shape)
{
    OCCGeometry* occgeo;
    occgeo = new OCCGeometry;
    BRep_Builder builder;
    TopoDS_Shape* s = (TopoDS_Shape*)shape;
    occgeo->shape = *s;

    occgeo->changed = 1;
    occgeo->BuildFMap();
    occgeo->CalcBoundingBox();
    //PrintContents(occgeo);
    return occgeo;
}

// Copy nodes and elements from origin_mesh into destination_mesh
void Cenos_MergeMesh(Ng_Mesh* orig_mesh, Ng_Mesh* dest_mesh, int index)
{
    Mesh* origin_mesh = (Mesh*)orig_mesh;
    Mesh* destination_mesh = (Mesh*)dest_mesh;

    //Geometries
    const OCCGeometry& origin_geom = *dynamic_pointer_cast<OCCGeometry>(origin_mesh->GetGeometry());
    const OCCGeometry& destination_geom = *dynamic_pointer_cast<OCCGeometry> (destination_mesh->GetGeometry());
    double eps = 1e-6 * destination_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 == 0)
            edge_index = destination_geom.emap.FindIndex(origin_geom.emap.FindKey(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 == 0)
            destination_mesh->AddFaceDescriptor(FaceDescriptor(destination_mesh->GetNFD() + 1, 1, 0, 1));
        else
            destination_mesh->AddFaceDescriptor(FaceDescriptor(index, 1, 0, 1));
        origin_mesh_dimension = 2;
    }
    else if (origin_geom.somap.Extent() == 1)
    {
        if (index == 0)
            solid_index = destination_geom.somap.FindIndex(origin_geom.somap.FindKey(1));
        else
            solid_index = index;
        origin_mesh_dimension = 3;
    }
    else // not supported
    {
        std::string message = "Ng_MergeMesh() does 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<PointIndex, PointIndex> node_ids;

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

        if (!exists)
        {
            destination_mesh->AddPoint(p);
            node_ids[opi] = PointIndex(destination_mesh->GetNP());
        }
    }

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

    //edge elements
    if (origin_mesh_dimension == 1)
    {
        for (int i = 1; i <= nEdgeEl; i++)
        {
            Segment seg = origin_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 = destination_mesh->GetNSeg() + 1;
            seg.cd2i = -1;

            destination_mesh->AddSegment(seg);
        }
    }

    //surface elements
    if (origin_mesh_dimension == 2)
    {
        for (int i = 1; i <= nFaceEl; i++)
        {
            //TODO check element type!!!
            Element2d el = origin_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() == QUAD)
                el.PNum(4) = node_ids[el.PNum(4)];

            el.SetIndex(destination_mesh->GetNFD());
            destination_mesh->AddSurfaceElement(el);
        }
    }

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

// Copy segments from origin_mesh into destination_mesh, assign edge index to segments
void Ng_CopySegments(Ng_Mesh* orig_mesh, Ng_Mesh* dest_mesh, int edge_index)
{
    Mesh* origin_mesh = (Mesh*)orig_mesh;
    Mesh* destination_mesh = (Mesh*)dest_mesh;

    //Geometries
    const OCCGeometry& origin_geom = *dynamic_pointer_cast<OCCGeometry>(origin_mesh->GetGeometry());
    const OCCGeometry& destination_geom = *dynamic_pointer_cast<OCCGeometry> (destination_mesh->GetGeometry());
    double eps = 1e-6 * destination_geom.GetBoundingBox().Diam();


    // 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<PointIndex, PointIndex> node_ids;

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

        if (!exists)
        {
            destination_mesh->AddPoint(p);
            node_ids[opi] = PointIndex(destination_mesh->GetNP());
        }
    }

    int nEdgeEl = origin_mesh->GetNSeg();

    //copy segments

    for (int i = 1; i <= nEdgeEl; i++)
    {
        Segment seg = origin_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 = destination_mesh->GetNSeg() + 1;
        seg.cd2i = -1;

        destination_mesh->AddSegment(seg);
    }
}

// Copy surface elements from origin_mesh into destination_mesh, assign edge index to segments
void Ng_CopySurfaceElements(Ng_Mesh* orig_mesh, Ng_Mesh* dest_mesh, int face_index)
{
    Mesh* origin_mesh = (Mesh*)orig_mesh;
    Mesh* destination_mesh = (Mesh*)dest_mesh;

    //Geometries
    const OCCGeometry& origin_geom = *dynamic_pointer_cast<OCCGeometry>(origin_mesh->GetGeometry());
    const OCCGeometry& destination_geom = *dynamic_pointer_cast<OCCGeometry> (destination_mesh->GetGeometry());
    double eps = 1e-6 * destination_geom.GetBoundingBox().Diam();


    // 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<PointIndex, PointIndex> node_ids;

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

    int nFaceEl = origin_mesh->GetNSE();

    // copy surface elements
    for (int i = 1; i <= nFaceEl; i++)
    {
        //TODO check element type!!!
        Element2d el = origin_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() == QUAD)
            el.PNum(4) = node_ids[el.PNum(4)];
        el.SetIndex(face_index);
        destination_mesh->AddSurfaceElement(el);
    }
}

// Prepare surface meshing
void Ng_PrepareSurfaceMeshing(Ng_Mesh* mesh)
{
    Mesh* occ_mesh = (Mesh*)mesh;

    const OCCGeometry& occgeom = *dynamic_pointer_cast<OCCGeometry> (occ_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 (SegmentIndex si = 0; si < occ_mesh->GetNSeg(); si++)
        {
            if ((*occ_mesh)[si].epgeominfo[1].edgenr == e)
            {
                gp_Pnt2d p2d;
                uvc++;
                p2d = cof->Value((*occ_mesh)[si].epgeominfo[0].dist);

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

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

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

        }

    }
}

// Reverse all segments in mesh
void Ng_ReverseSegments(Ng_Mesh* mesh)
{
    Mesh* occ_mesh = (Mesh*)mesh;

    for (SegmentIndex si = 0; si < occ_mesh->GetNSeg(); si++)
    {
        swap((*occ_mesh)[si][0], (*occ_mesh)[si][1]);
        swap((*occ_mesh)[si].epgeominfo[0].dist, (*occ_mesh)[si].epgeominfo[1].dist);
        swap((*occ_mesh)[si].epgeominfo[0].edgenr, (*occ_mesh)[si].epgeominfo[1].edgenr);
        swap((*occ_mesh)[si].epgeominfo[0].u, (*occ_mesh)[si].epgeominfo[1].u);
        swap((*occ_mesh)[si].epgeominfo[0].v, (*occ_mesh)[si].epgeominfo[1].v);
    }
}

// Reverse all faces in mesh
void Ng_ReverseFaces(Ng_Mesh* mesh)
{
    Mesh* occ_mesh = (Mesh*)mesh;

    for (SurfaceElementIndex si = 0; si < occ_mesh->GetNSE(); si++)
    {
        (*occ_mesh)[si].Invert();
    }
}

Ng_LocalH Ng_GetLocalH(Ng_Mesh* mesh)
{
    ((Mesh*)mesh)->CalcLocalHFromPointDistances(0.4);
    return (Ng_LocalH)((Mesh*)mesh)->GetLocalH().get();
}

void Ng_CopyLocalH(Ng_Mesh* orig_mesh, Ng_Mesh* dest_mesh)
{
    Mesh* origin_mesh = (Mesh*)orig_mesh;
    Mesh* destination_mesh = (Mesh*)dest_mesh;
    destination_mesh->SetLocalH(origin_mesh->GetLocalH());
}

double Ng_GetMaxH(Ng_Mesh* occ_mesh)
{
    Mesh* mesh = (Mesh*)occ_mesh;
    mesh->CalcLocalHFromPointDistances(0.4);

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

void Ng_ExportMeshToGmesh2(Ng_Mesh* mesh, const char* filename)
{
    //Mesh* m = (Mesh*)mesh;
    Cenos_WriteGmsh2Format(*(Mesh*)mesh, string(filename));
}

void Ng_ExportMeshToOpenFOAM(Ng_Mesh* mesh, const char* filename)
{
    //Mesh* m = (Mesh*)mesh;
    WriteOpenFOAM15xFormat(*(Mesh*)mesh, string(filename), false);
}

// Extract the solid map from the OCC geometry
// The solid map basically gives an index to each solid in the geometry, 
// which can be used to access a specific solid
Ng_Result Ng_GetSolidMap(Ng_OCC_Geometry* geom,
    Ng_OCC_TopTools_IndexedMapOfShape* SolidMap)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    TopTools_IndexedMapOfShape* occsomap = (TopTools_IndexedMapOfShape*)SolidMap;

    // Copy the solid map from the geometry to the given variable
    occsomap->Assign(occgeom->somap);

    if (occsomap->Extent())
    {
        return NG_OK;
    }
    else
    {
        return NG_ERROR;
    }
}

// Extract the face map from the OCC geometry
// The face map basically gives an index to each face in the geometry, 
// which can be used to access a specific face
// copy of Ng_OCC_GetFMap
Ng_Result Ng_GetFaceMap(Ng_OCC_Geometry* geom,
    Ng_OCC_TopTools_IndexedMapOfShape* FaceMap)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    TopTools_IndexedMapOfShape* occfmap = (TopTools_IndexedMapOfShape*)FaceMap;

    // Copy the face map from the geometry to the given variable
    occfmap->Assign(occgeom->fmap);

    if (occfmap->Extent())
    {
        return NG_OK;
    }
    else
    {
        return NG_ERROR;
    }
}

// Extract the edge map from the OCC geometry
// The edge map basically gives an index to each edge in the geometry, 
// which can be used to access a specific edge
Ng_Result Ng_GetEdgeMap(Ng_OCC_Geometry* geom,
    Ng_OCC_TopTools_IndexedMapOfShape* EdgeMap)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    TopTools_IndexedMapOfShape* occedgemap = (TopTools_IndexedMapOfShape*)EdgeMap;

    // Copy the edge map from the geometry to the given variable
    occedgemap->Assign(occgeom->emap);

    if (occedgemap->Extent())
    {
        return NG_OK;
    }
    else
    {
        return NG_ERROR;
    }
}

int Ng_GetEdgeElementIndex(Ng_Mesh* mesh, int num)
{
    //        if (((Mesh*)mesh)->GetDimension() == 3)
    //            return ((Mesh*)mesh)->LineSegment(num).edgenr;
    //        else
    return ((Mesh*)mesh)->LineSegment(num).epgeominfo[0].edgenr;
}

int Ng_GetSurfaceElementIndex(Ng_Mesh* mesh, int num)
{
    int fd_id = ((Mesh*)mesh)->SurfaceElement(num).GetIndex();
    return ((Mesh*)mesh)->GetFaceDescriptor(fd_id).SurfNr();
}

int Ng_GetVolumeElementIndex(Ng_Mesh* mesh, int num)
{
    return ((Mesh*)mesh)->VolumeElement(num).GetIndex();
}

void Ng_RedirectCout(void* ptr_filestream)
{
    mycout = (ofstream*)ptr_filestream;
}

void Ng_CalculateSurfacesOfNode(Ng_Mesh* mesh)
{
    Mesh* m = (Mesh*)mesh;
    m->CalcSurfacesOfNode();
}

void Ng_DumpSegments(Ng_Mesh* mesh, void* ptr_filestream)
{
    Mesh* m = (Mesh*)mesh;
    ostream* file = (ofstream*)ptr_filestream;
    (*file) << "GetNCD2Names() GetNFD() " << std::endl;

    (*file) << m->GetNCD2Names() << " " << m->GetNFD() << std::endl;

    (*file) << "pp1 p2 edgenr singedge_left singedge_right trignum[0] u0 v0 trignum[1] u1 v1 dist1 dist2 si cd2i domin domout tlosurf surfnr1 surfnr2" << std::endl;

    for (int i = 1; i <= m->GetNSeg(); i++)
    {
        Segment seg = m->LineSegment(i);
        (*file) << seg[0] << " " << seg[1] << " " << seg.edgenr << " " << seg.singedge_left <<
            " " << seg.singedge_right << " " << seg.geominfo[0].trignum << " " << seg.epgeominfo[0].u << " " << seg.epgeominfo[0].v << " " <<
            seg.geominfo[1].trignum << " " << seg.epgeominfo[1].u << " " << seg.epgeominfo[1].v << " " << seg.epgeominfo[0].dist << " " << seg.epgeominfo[1].dist
            << " " << seg.si << " " << seg.cd2i << " " << seg.domin << " " << seg.domout << " " << seg.tlosurf << " " <<
            seg.surfnr1 << " " << seg.surfnr2 << std::endl;
    }
}

int Ng_GetNSeg(Ng_Mesh* mesh)
{
    return ((Mesh*)mesh)->GetNSeg();
}

int Ng_GetSegment(Ng_Mesh* mesh, int num, int* pi, int* matnum)
{
    const Segment& seg = ((Mesh*)mesh)->LineSegment(num);
    pi[0] = seg[0];
    pi[1] = seg[1];

    if (matnum)
        *matnum = seg.edgenr;
    if (seg[2] == 0)
    {
        return 1;
    }
    else
    {
        pi[2] = seg[2];
        return 8;
    }
}

// Manually add a segment element of a given type to an existing mesh object
void Ng_AddSegmentElement(Ng_Mesh* mesh, int pi1, int pi2, int edgeIndex)
{
    Mesh* m = (Mesh*)mesh;
    const Point3d& p1 = m->Point(pi1);
    const Point3d& p2 = m->Point(pi2);

    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 = 1;
    seg.edgenr = m->GetNSeg() + 1;

    m->AddSegment(seg);
}

// Set a local limit on the maximum mesh size by copyin max size from origin mesh
void Ng_RestrictMeshSizeMesh(Ng_Mesh* orig_mesh, Ng_Mesh* dest_mesh)
{
    Mesh* origin_mesh = (Mesh*)orig_mesh;
    Mesh* destination_mesh = (Mesh*)dest_mesh;

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

// Set a local limit on the maximum mesh size by copyin max size from origin mesh
void Ng_RestrictMeshSizeLocalH(Ng_Mesh* occ_mesh, Ng_LocalH* localh)
{
    Mesh* mesh = (Mesh*)occ_mesh;
    for (PointIndex opi : mesh->Points().Range())
    {
        const Point3d& p = (*mesh)[opi];
        mesh->RestrictLocalH(p, ((LocalH*)localh)->GetH(p));
    }
}

Ng_Result Ng_GenerateBoundaryLayer(Ng_Mesh* mesh,
    int* surfid_arr, int surfid_count,
    double* heights_arr, int heights_count)
{
    Mesh* m = (Mesh*)mesh;

    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;

    try
    {
        GenerateBoundaryLayer(*m, blp);
        //   RemoveIllegalElements(*m);
        OptimizeVolume(mparam, *m);
    }
    catch (NgException e)
    {
        std::cout << "Netgen Exception: " << e.what() << std::endl;
        return NG_ERROR;
    }

    return NG_OK;
}

Ng_Result Ng_GenerateBoundaryLayer2(Ng_Mesh* mesh,
    int dom_nr, double* heights_arr, int heights_count)
{
    Mesh* m = (Mesh*)mesh;

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

    try
    {
        GenerateBoundaryLayer2(*m, dom_nr, thicknesses, false);
    }
    catch (NgException e)
    {
        std::cout << "Netgen Exception: " << e.what() << std::endl;
        return NG_ERROR;
    }

    return NG_OK;
}

void Ng_AddFaceDescriptor(Ng_Mesh* mesh, int faceId, int dominId, int domOutId, int TLOface)
{
    Mesh* m = (Mesh*)mesh;
    m->AddFaceDescriptor(FaceDescriptor(faceId, dominId, domOutId, TLOface));
}

void Ng_AddEdgeDescriptor(Ng_Mesh* mesh, int edgeId)
{
    Mesh* m = (Mesh*)mesh;
    EdgeDescriptor  ed;
    ed.SetTLOSurface(edgeId);
    m->AddEdgeDescriptor(ed);
}

// Manually add a surface element of a given type to an existing mesh object with
// added surface index
void Ng_AddSurfaceElement(Ng_Mesh* mesh, Ng_Surface_Element_Type et,
    int* pi, int surfIndx)
{
    Mesh* m = (Mesh*)mesh;
    Element2d el(3);
    if (et == NG_TRIG)
        el = Element2d(3);
    else if (et == NG_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 == NG_QUAD)
        el.PNum(4) = pi[3];

    m->AddSurfaceElement(el);
}

// Manually add a volume element of a given type to an existing mesh object
// added volume index
// added prism check
void Ng_AddVolumeElement(Ng_Mesh* mesh, Ng_Volume_Element_Type et,
    int* pi, int volIndx)
{
    Mesh* m = (Mesh*)mesh;
    int nodeCount = 4;
    if (et == Ng_Volume_Element_Type::NG_PRISM)
    {
        nodeCount = 6;
    }

    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 == Ng_Volume_Element_Type::NG_PRISM)
    {
        el.PNum(5) = pi[4];
        el.PNum(6) = pi[5];
    }

    m->AddVolumeElement(el);
}

// --------------------- OCC Geometry / Meshing Utility Functions -------------------

Ng_OCC_Geometry* Ng_OCC_ShapeToGeometry(Ng_TopoDS_Shape* s)
{
    OCCGeometry* geom = CreateOCCGeometryFromTopoDS((Cenos_TopoDS_Shape*)s);
    return (Ng_OCC_Geometry*)geom;
}

// Locally limit the size of the mesh to be generated at various points 
// based on the topology of the geometry
// copy of Ng_OCC_SetLocalMeshSize
Ng_Result Ng_SetLocalMeshSize(Ng_OCC_Geometry* geom,
    Ng_Mesh* mesh,
    Ng_Meshing_Parameters* mp)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    Mesh* me = (Mesh*)mesh;
    me->SetGeometry(shared_ptr<OCCGeometry>(occgeom, &NOOP_Deleter));

    me->geomtype = Mesh::GEOM_OCC;

    mp->Transfer_Parameters();

    if (mp->closeedgeenable)
        mparam.closeedgefac = mp->closeedgefact;

    // Delete the mesh structures in order to start with a clean 
    // slate
    me->DeleteMesh();

    OCCSetLocalMeshSize(*occgeom, *me, mparam, occparam);

    return(NG_OK);
}

// Set the geometry / topology of mesh
void Ng_SetGeometry(Ng_OCC_Geometry* geom, Ng_Mesh* mesh)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    Mesh* me = (Mesh*)mesh;
    me->SetGeometry(shared_ptr<OCCGeometry>(occgeom, &NOOP_Deleter));
}

// Mesh single edgeId
Ng_Result Ng_DivideEdges(Ng_OCC_Geometry* geom,
    Ng_Mesh* mesh,
    Ng_Meshing_Parameters* meshparams)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    Mesh* me = (Mesh*)mesh;
    me->SetGeometry(shared_ptr<OCCGeometry>(occgeom, &NOOP_Deleter));
    meshparams->Transfer_Parameters();
    double eps = 1e-6 * occgeom->GetBoundingBox().Diam();

    // add vertices to mesh as points
    for (int i = 1; i <= occgeom->vmap.Extent(); i++)
    {
        gp_Pnt pnt = BRep_Tool::Pnt(TopoDS::Vertex(occgeom->vmap(i)));
        MeshPoint mpnt(Point<3>(pnt.X(), pnt.Y(), pnt.Z()));

        bool exists = 0;
        for (PointIndex pi : me->Points().Range())
            if (Dist2((*me)[pi], Point<3>(mpnt)) < eps * eps)
            {
                exists = true;
                break;
            }

        if (!exists)
            me->AddPoint(mpnt);
    }


    PointIndex first_ep = me->Points().Range().end();
    auto vertexrange = me->Points().Range();


    for (int edge_i = 1; edge_i <= occgeom->emap.Extent(); edge_i++)
    {
        TopoDS_Edge edge = TopoDS::Edge(occgeom->emap(edge_i));
        if (BRep_Tool::Degenerated(edge))
        {
            //(*testout) << "ignoring degenerated edge" << endl;
            continue;
        }

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

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

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


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

        Cenos_DivideEdge(edge, mp, params, *me, mparam);


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

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

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

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

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

        (*testout) << "NP = " << me->GetNP() << endl;

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

            seg[0] = pnums[i - 1];
            seg[1] = pnums[i];
            seg.edgenr = me->GetNSeg() + 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);
            }

            me->AddSegment(seg);
        }
    }


    if ((me->GetNP()) && (me->GetNSeg()))
    {
        return NG_OK;
    }
    else
    {
        return NG_ERROR;
    }
}

// Mesh the edges and add Face descriptors to prepare for surface meshing
// copy of Ng_OCC_GenerateEdgeMesh
Ng_Result Ng_GenerateEdgeMesh(Ng_OCC_Geometry* geom,
    Ng_Mesh* mesh,
    Ng_Meshing_Parameters* mp)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    Mesh* me = (Mesh*)mesh;
    me->SetGeometry(shared_ptr<OCCGeometry>(occgeom, &NOOP_Deleter));

    mp->Transfer_Parameters();

    OCCFindEdges(*occgeom, *me, mparam);

    if ((me->GetNP()) && (me->GetNFD()))
    {
        return NG_OK;
    }
    else
    {
        return NG_ERROR;
    }
}

// Mesh the edges and add Face descriptors to prepare for surface meshing
// copy of Ng_OCC_GenerateSurfaceMesh
Ng_Result Ng_GenerateSurfaceMesh(Ng_OCC_Geometry* geom,
    Ng_Mesh* mesh,
    Ng_Meshing_Parameters* mp)
{
    int numpoints = 0;
    (*mycout) << "Starting GenerateSurfaceMesh" << endl << flush;

    OCCGeometry* occgeom = (OCCGeometry*)geom;
    Mesh* me = (Mesh*)mesh;
    if (!me->GetNFD())
        me->AddFaceDescriptor(FaceDescriptor(1, 1, 0, 0));

    me->SetGeometry(shared_ptr<OCCGeometry>(occgeom, &NOOP_Deleter));
    (*mycout) << "Geometry is set" << endl << flush;

    // Set the internal meshing parameters structure from the nglib meshing 
    // parameters structure
    mp->Transfer_Parameters();
    (*mycout) << "Parameters are transferred" << endl << flush;

    // Only go into surface meshing if the face descriptors have already been added
    if (!me->GetNFD())
        return NG_ERROR;

    numpoints = me->GetNP();

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

    // Check and if required, enable surface mesh optimisation step
    if (mp->optsurfmeshenable)
    {
        perfstepsend = MESHCONST_OPTSURFACE;
    }
    (*mycout) << "Start mesh surface" << endl << flush;
    try
    {
        OCCMeshSurface(*occgeom, *me, mparam);

        OCCOptimizeSurface(*occgeom, *me, mparam);

    }
    catch (NgException e)
    {
        std::cout << "Netgen Exception: " << e.what() << std::endl;
        return NG_ERROR;
    }
    (*mycout) << "Done meshing surface" << endl << flush;
    (*mycout) << "Optimizing surface" << endl << flush;

    (*mycout) << "Done optimizingsurface" << endl << flush;

    me->CalcSurfacesOfNode();
    (*mycout) << "Done CalcSurfacesOfNode" << endl << flush;

    (*mycout) << me->GetNP() << endl << flush;
    (*mycout) << numpoints << endl << flush;
    (*mycout) << me->GetNSE() << endl << flush;

    if (mp->second_order)
        me->GetGeometry()->GetRefinement().MakeSecondOrder(*me);

    if (me->GetNP() < numpoints)
        return NG_ERROR;

    if (me->GetNSE() <= 0)
        return NG_ERROR;

    return NG_OK;
}

// copy of Ng_OCC_Uniform_Refinement
void Ng_Refine(Ng_Mesh* mesh)
{

    /*Mesh* me = (Mesh*)mesh;
    me->GetGeometry()->GetRefinement().Refine(*me);
    me->UpdateTopology();*/


    BisectionOptions biopt;
    biopt.usemarkedelements = 1;
    biopt.refine_p = 0;
    biopt.refine_hp = 0;
    Mesh* me = (Mesh*)mesh;

    for (int i = 1; i <= me->GetNSE(); i++)
        me->SurfaceElement(i).SetRefinementFlag(false);

    me->GetGeometry()->GetRefinement().Bisect(*me, biopt);
    //me->GetGeometry()->GetRefinement().Refine(*me);
    me->UpdateTopology();
    me->GetCurvedElements().SetIsHighOrder(false);

    /* other possible option - check if same result!
    * ( (OCCGeometry*)geom ) -> GetRefinement().Refine ( * (Mesh*) mesh );
   */
}

// Set refinement flag for element
void Ng_SetElementRefinement(Ng_Mesh* mesh, int el_index, bool flag)
{
    Mesh* me = (Mesh*)mesh;
    if (me->GetDimension() == 3)
    {
        me->VolumeElement(el_index).SetRefinementFlag(flag);
    }
    else
    {
        me->SurfaceElement(el_index).SetRefinementFlag(flag);
    }
}

void Cenos_ExportMeshToGmesh2(Ng_Mesh* mesh, const char* filename)
{
    //Mesh* m = (Mesh*)mesh;
    Cenos_WriteGmsh2Format(*(Mesh*)mesh, string(filename));
    //netgen::WriteUserFormat(string("Cenos Gmsh2 Format"), *(Mesh*)mesh, /* geom, */ string(filename));
}

void Cenos_GenerateBoundaryLayer(Ng_Mesh* mesh,
    int* surfid_arr, int surfid_count,
    double* heights_arr, int heights_count)
{
    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*)mesh, blp);
    std::cout << "Meshing BL done" << std::endl;

}

Ng_OCC_Geometry* Cenos_OCC_ShapeToGeometry(Ng_TopoDS_Shape* s)
{
    OCCGeometry* geom = CreateOCCGeometryFromTopoDS((Cenos_TopoDS_Shape*)s);
    return (Ng_OCC_Geometry*)geom;
}

Ng_Result Cenos_OCC_GetSoMap(Ng_OCC_Geometry* geom,
    Ng_OCC_TopTools_IndexedMapOfShape* SoMap)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    TopTools_IndexedMapOfShape* occsomap = (TopTools_IndexedMapOfShape*)SoMap;

    // Copy the face map from the geometry to the given variable
    occsomap->Assign(occgeom->somap);

    if (occsomap->Extent())
    {
        return NG_OK;
    }
    else
    {
        return NG_ERROR;
    }
}

Ng_Result Cenos_OCC_GetEdgeMap(Ng_OCC_Geometry* geom,
    Ng_OCC_TopTools_IndexedMapOfShape* EdgeMap)
{
    OCCGeometry* occgeom = (OCCGeometry*)geom;
    TopTools_IndexedMapOfShape* occedgemap = (TopTools_IndexedMapOfShape*)EdgeMap;

    // Copy the face map from the geometry to the given variable
    occedgemap->Assign(occgeom->emap);

    if (occedgemap->Extent())
    {
        return NG_OK;
    }
    else
    {
        return NG_ERROR;
    }
}

int Cenos_GetEdgeElementIndex(Ng_Mesh* mesh, int num)
{
    //        if (((Mesh*)mesh)->GetDimension() == 3)
    //            return ((Mesh*)mesh)->LineSegment(num).edgenr;
    //        else
    return ((Mesh*)mesh)->LineSegment(num).epgeominfo[0].edgenr;
}

int Cenos_GetSurfaceElementIndex(Ng_Mesh* mesh, int num)
{
    return ((Mesh*)mesh)->SurfaceElement(num).GetIndex();
}


int Cenos_GetVolumeElementIndex(Ng_Mesh* mesh, int num)
{
    return ((Mesh*)mesh)->VolumeElement(num).GetIndex();
}

void Cenos_RedirectCout(void* ptr_filestream)
{
    mycout = (ofstream*)ptr_filestream;
}

void Cenos_DumpSegments(Ng_Mesh* mesh, void* ptr_filestream)
{
    Mesh* m = (Mesh*)mesh;
    ostream* file = (ofstream*)ptr_filestream;

    (*file) << m->GetNCD2Names() << " " << m->GetNFD() << std::endl;

    m->CalcSurfacesOfNode();
    (*file) << m->GetNOpenSegments() << " " << m->GetNOpenElements() << " " << m->CheckConsistentBoundary() << " " << m->CheckOverlappingBoundary() << std::endl;
    int nrOfSeg = m->GetNSeg();
    for (int i = 1; i <= nrOfSeg; i++)
    {
        Segment seg = m->LineSegment(i);
        (*file) << seg[0] << " " << seg[1] << " " << seg.edgenr << " " << seg.singedge_left << " " << seg.singedge_right << " " << seg.geominfo[0].trignum << " " << seg.geominfo[0].u << " " << seg.geominfo[0].v << seg.geominfo[1].trignum << " " << seg.geominfo[1].u << " " << seg.geominfo[1].v << std::endl;


    }
}

// Manually add a segment element of a given type to an existing mesh object
void Cenos_AddSegmentElement(Ng_Mesh* mesh, int pi1, int pi2, int edgeIndex, double* zeroNode)
{
    Mesh* m = (Mesh*)mesh;
    const Point3d& p1 = m->Point(pi1);
    const Point3d& p2 = m->Point(pi2);

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

    seg.epgeominfo[0].dist = sqrt(pow(p1.X() - zeroNode[0], 2) +
        pow(p1.Y() - zeroNode[1], 2) +
        pow(p1.Z() - zeroNode[2], 2));
    seg.epgeominfo[1].dist = sqrt(pow(p2.X() - zeroNode[0], 2) +
        pow(p2.Y() - zeroNode[1], 2) +
        pow(p2.Z() - zeroNode[2], 2));
    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 = 1;
    seg.edgenr = m->GetNSeg() + 1;

    m->AddSegment(seg);
}

// Manually add a surface element of a given type to an existing mesh object
void Cenos_AddSurfaceElementUV(Ng_Mesh* mesh, Ng_Surface_Element_Type et,
    int* pi, int surfIndx, double* uv1, double* uv2, double* uv3)
{
    Mesh* m = (Mesh*)mesh;
    Element2d el(3);
    el.SetIndex(surfIndx);
    el.PNum(1) = pi[0];
    el.PNum(2) = pi[1];
    el.PNum(3) = pi[2];
    el.GeomInfoPi(1).u = uv1[0];
    el.GeomInfoPi(1).v = uv1[1];
    el.GeomInfoPi(2).u = uv2[0];
    el.GeomInfoPi(2).v = uv2[1];
    el.GeomInfoPi(3).u = uv3[0];
    el.GeomInfoPi(3).v = uv3[1];
    m->AddSurfaceElement(el);
}

// Manually add a surface element of a given type to an existing mesh object
void Cenos_AddSurfaceElement(Ng_Mesh* mesh, Ng_Surface_Element_Type et,
    int* pi, int surfIndx)
{
    Mesh* m = (Mesh*)mesh;
    Element2d el(3);
    el.SetIndex(surfIndx);
    el.PNum(1) = pi[0];
    el.PNum(2) = pi[1];
    el.PNum(3) = pi[2];
    m->AddSurfaceElement(el);
}

// Manually add a volume element of a given type to an existing mesh object
void Cenos_AddVolumeElement(Ng_Mesh* mesh, Ng_Volume_Element_Type et,
    int* pi, int volIndx)
{
    Mesh* m = (Mesh*)mesh;
    int nodeCount = 4;
    if (et == Ng_Volume_Element_Type::NG_PRISM)
    {
        nodeCount = 6;
    }

    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 == Ng_Volume_Element_Type::NG_PRISM)
    {
        el.PNum(5) = pi[4];
        el.PNum(6) = pi[5];
    }

    m->AddVolumeElement(el);
}

Ng_Mesh* Cenos_NewMesh()
{
    Mesh* mesh = new Mesh;
    return (Ng_Mesh*)(void*)mesh;
}

void Cenos_AddFaceDescriptor(Ng_Mesh* mesh, int faceId, int dominId, int domOutId, int TLOface)
{
    Mesh* m = (Mesh*)mesh;
    m->AddFaceDescriptor(FaceDescriptor(faceId, dominId, domOutId, TLOface));
}


void Cenos_AddEdgeDescriptor(Ng_Mesh* mesh, int edgeId)
{
    Mesh* m = (Mesh*)mesh;
    EdgeDescriptor  ed;
    ed.SetTLOSurface(edgeId);
    m->AddEdgeDescriptor(ed);
}

/*! \file writegmsh2.cpp
*  \brief Export Netgen Mesh in the GMSH v2.xx File format
*  \author Philippose Rajan
*  \date 02 November 2008
*
*  This function extends the export capabilities of
*  Netgen to include the GMSH v2.xx File Format.
*
*  Current features of this function include:
*
*  1. Exports Triangles, Quadrangles and Tetrahedra \n
*  2. Supports upto second order elements of each type
*
*/


/*! GMSH v2.xx mesh format export function
*
*  This function extends the export capabilities of
*  Netgen to include the GMSH v2.xx File Format.
*
*  Current features of this function include:
*
*  1. Exports Triangles, Quadrangles and Tetrahedra \n
*  2. Supports upto second order elements of each type
*
*/
void Cenos_WriteGmsh2Format(const Mesh& mesh,
    const string& filename)
{
    ofstream outfile(filename.c_str());
    outfile.precision(6);
    outfile.setf(ios::fixed, ios::floatfield);
    outfile.setf(ios::showpoint);

    int np = mesh.GetNP();  /// number of points in mesh
    int ne = mesh.GetNE();  /// number of 3D elements in mesh
    int nse = mesh.GetNSE();  /// number of surface elements (BC)
    int i, j, k, l;


    /*
    * 3D section : Volume elements (currently only tetrahedra)
    */

    if ((ne > 0)
        && (mesh.VolumeElement(1).GetNP() <= 10)
        && (mesh.SurfaceElement(1).GetNP() <= 6))
    {
        cout << "Write GMSH v2.xx Format \n";
        cout << "The GMSH v2.xx export is currently available for elements upto 2nd Order\n" << endl;

        int inverttets = mparam.inverttets;
        int invertsurf = mparam.inverttrigs;

        /// Prepare GMSH 2.0 file (See GMSH 2.0 Documentation)
        outfile << "$MeshFormat\n";
        outfile << (float)2.0 << " "
            << (int)0 << " "
            << (int)sizeof(double) << "\n";
        outfile << "$EndMeshFormat\n";


        // write physical entities
        int maxSurfIndex = 0;
        for (i = 1; i <= nse; i++)
        {
            int elType = 0;
            Element2d el = mesh.SurfaceElement(i);
            if (el.GetType() == ELEMENT_TYPE::TRIG) elType = GMSH_TRIG;	//// GMSH Type for a 3 node triangle
            if (el.GetType() == ELEMENT_TYPE::QUAD) elType = GMSH_QUAD;	//// GMSH Type for a 6 node prism
            if (el.GetType() == ELEMENT_TYPE::TRIG6) elType = GMSH_TRIG6;  //// GMSH Type for a 6 node triangle
            if (elType == 0)
            {
                cout << " Invalid surface element type for Gmsh 2.0 3D-Mesh Export Format !\n";
                return;
            }

            int elIndex = mesh.GetFaceDescriptor(el.GetIndex()).BCProperty();
            if (maxSurfIndex < elIndex)
                maxSurfIndex = elIndex;
        }

        int maxSolidIndex = 0;
        for (i = 1; i <= ne; i++)
        {
            int elType = 0;

            Element el = mesh.VolumeElement(i);
            if (inverttets) el.Invert();

            if (el.GetType() == ELEMENT_TYPE::TET) elType = GMSH_TET;    //// GMSH Element type for 4 node tetrahedron
            if (el.GetType() == ELEMENT_TYPE::PRISM) elType = GMSH_PRISM;    //// GMSH Element type for 6 node prism
            if (el.GetType() == ELEMENT_TYPE::TET10) elType = GMSH_TET10;    //// GMSH Element type for 10 node tetrahedron
            if (elType == 0)
            {
                cout << " Invalid volume element type for Gmsh 2.0 3D-Mesh Export Format !\n";
                return;
            }
            int elIndex = el.GetIndex();
            if (maxSolidIndex < elIndex)
                maxSolidIndex = elIndex;
        }

        /// Write nodes
        outfile << "$Nodes\n";
        outfile << np << "\n";

        for (i = 1; i <= np; i++)
        {
            const Point3d& p = mesh.Point(i);
            outfile << i << " "; /// node number
            outfile << p.X() << " ";
            outfile << p.Y() << " ";
            outfile << p.Z() << "\n";
        }

        outfile << "$EndNodes\n";

        /// write elements (both, surface elements and volume elements)
        outfile << "$Elements\n";
        outfile << ne + nse << "\n";  ////  number of elements + number of surfaces BC

        for (i = 1; i <= nse; i++)
        {
            int elType = 0;

            Element2d el = mesh.SurfaceElement(i);
            if (invertsurf) el.Invert();

            if (el.GetType() == ELEMENT_TYPE::TRIG) elType = GMSH_TRIG;	//// GMSH Type for a 3 node triangle
            if (el.GetType() == ELEMENT_TYPE::QUAD) elType = GMSH_QUAD;	//// GMSH Type for a 6 node prism
            if (el.GetType() == ELEMENT_TYPE::TRIG6) elType = GMSH_TRIG6;  //// GMSH Type for a 6 node triangle
            if (elType == 0)
            {
                cout << " Invalid surface element type for Gmsh 2.0 3D-Mesh Export Format !\n";
                return;
            }

            outfile << i;
            outfile << " ";
            outfile << elType;
            outfile << " ";
            outfile << "2";                  //// Number of tags (2 => Physical and elementary entities)
            outfile << " ";
            outfile << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() << " ";
            /// that means that physical entity = elementary entity (arbitrary approach)
            outfile << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() << " ";
            for (j = 1; j <= el.GetNP(); j++)
            {
                outfile << " ";
                if ((elType == GMSH_TRIG) || (elType == GMSH_TRIG6))
                {
                    outfile << el.PNum(triGmsh[l]);
                }
                else if ((elType == GMSH_QUAD) || (elType == GMSH_QUAD8))
                {
                    outfile << el.PNum(quadGmsh[l]);
                }
            }
            outfile << "\n";
        }


        for (i = 1; i <= ne; i++)
        {
            int elType = 0;

            Element el = mesh.VolumeElement(i);
            if (inverttets) el.Invert();

            if (el.GetType() == ELEMENT_TYPE::TET) elType = GMSH_TET;    //// GMSH Element type for 4 node tetrahedron
            if (el.GetType() == ELEMENT_TYPE::PRISM) elType = GMSH_PRISM;    //// GMSH Element type for 4 node tetrahedron
            if (el.GetType() == ELEMENT_TYPE::TET10) elType = GMSH_TET10;    //// GMSH Element type for 4 node tetrahedron
            if (elType == 0)
            {
                cout << " Invalid volume element type for Gmsh 2.0 3D-Mesh Export Format !\n";
                return;
            }

            outfile << nse + i;                       //// element number (Remember to add on surface elements)
            outfile << " ";
            outfile << elType;
            outfile << " ";
            outfile << "2";                   //// Number of tags (2 => Physical and elementary entities)
            outfile << " ";
            outfile << maxSurfIndex + el.GetIndex();
            /// that means that physical entity = elementary entity (arbitrary approach)
            outfile << " ";
            outfile << maxSurfIndex + el.GetIndex();   /// volume number
            outfile << " ";
            for (j = 1; j <= el.GetNP(); j++)
            {
                outfile << " ";
                if ((elType == GMSH_TET) || (elType == GMSH_TET10))
                {
                    outfile << el.PNum(tetGmsh[j]);
                }
                else if ((elType == GMSH_PRISM))
                {
                    outfile << el.PNum(prismGmsh[j]);
                }
            }
            outfile << "\n";
        }
        outfile << "$EndElements\n";
    }
    /*
    * End of 3D section
    */


    /*
    * 2D section : available for triangles and quadrangles
    *              upto 2nd Order
    */
    else if (ne == 0)   /// means that there's no 3D element
    {
        cout << "\n Write Gmsh v2.xx Surface Mesh (triangle and/or quadrangles upto 2nd Order)" << endl;

        /// Prepare GMSH 2.0 file (See GMSH 2.0 Documentation)
        outfile << "$MeshFormat\n";
        outfile << (float)2.0 << " "
            << (int)0 << " "
            << (int)sizeof(double) << "\n";
        outfile << "$EndMeshFormat\n";

        /// Write nodes
        outfile << "$Nodes\n";
        outfile << np << "\n";

        for (i = 1; i <= np; i++)
        {
            const Point3d& p = mesh.Point(i);
            outfile << i << " "; /// node number
            outfile << p.X() << " ";
            outfile << p.Y() << " ";
            outfile << p.Z() << "\n";
        }
        outfile << "$EndNodes\n";

        /// write triangles & quadrangles
        outfile << "$Elements\n";
        outfile << nse << "\n";

        for (k = 1; k <= nse; k++)
        {
            int elType = 0;

            const Element2d& el = mesh.SurfaceElement(k);

            if (el.GetNP() == 3) elType = GMSH_TRIG;   //// GMSH Type for a 3 node triangle
            if (el.GetNP() == 6) elType = GMSH_TRIG6;  //// GMSH Type for a 6 node triangle
            if (el.GetNP() == 4) elType = GMSH_QUAD;   //// GMSH Type for a 4 node quadrangle
            if (el.GetNP() == 8) elType = GMSH_QUAD8;  //// GMSH Type for an 8 node quadrangle
            if (elType == 0)
            {
                cout << " Invalid surface element type for Gmsh 2.0 2D-Mesh Export Format !\n";
                return;
            }

            outfile << k;
            outfile << " ";
            outfile << elType;
            outfile << " ";
            outfile << "2";
            outfile << " ";
            outfile << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() << " ";
            /// that means that physical entity = elementary entity (arbitrary approach)
            outfile << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() << " ";
            for (l = 1; l <= el.GetNP(); l++)
            {
                outfile << " ";
                if ((elType == GMSH_TRIG) || (elType == GMSH_TRIG6))
                {
                    outfile << el.PNum(triGmsh[l]);
                }
                else if ((elType == GMSH_QUAD) || (elType == GMSH_QUAD8))
                {
                    outfile << el.PNum(quadGmsh[l]);
                }
            }
            outfile << "\n";
        }
        outfile << "$EndElements\n";
    }
    /*
    * End of 2D section
    */

    else
    {
        cout << " Invalid element type for Gmsh v2.xx Export Format !\n";
    }
} // End: WriteGmsh2Format

void Cenos_DivideEdge(TopoDS_Edge& edge, netgen::NgArray<netgen::MeshPoint>& ps, netgen::NgArray<double>& params, Mesh& mesh, const netgen::MeshingParameters& mparam)
{
    double s0, s1;
    int nsubedges = 1;
    gp_Pnt pnt, oldpnt;
    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 (int i = 1; i <= tmpVal; i++)
    {
        oldpnt = pnt;
        pnt = c->Value(s0 + (i / double(DIVIDEEDGESECTIONS)) * (s1 - s0));
        hvalue[i] = hvalue[i - 1] +
            1.0 / mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) *
            pnt.Distance(oldpnt);

        //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))
        //	   <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;
    }

    //  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] = s0 + (i1 / double(DIVIDEEDGESECTIONS)) * (s1 - s0);
            pnt = c->Value(params[i]);
            ps[i - 1] = MeshPoint(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] = s0;
    params[nsubedges] = s1;

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