OCC + nglib: unv mesh with boundaries

Hello Victor,
we are using Netgen at our company, here is some code snippet of how we use it to produce a mesh from a TopoDS_Shape, it might not work as is but should be pretty close to working (I extract this code from a C++ class we use):

// setting of meshing parameters, tolerance,...
netgen::mparam.maxh = 1.; // some max edge length
netgen::mparam.minh = 1.; // some min edge length
netgen::mparam.optimize2d = "smsmsmSmSmSm";
netgen::mparam.grading = 0.3;
netgen::mparam.giveuptol2d=3000;
netgen::mparam.optsteps2d=3;

nglib::Ng_Init();

// -------------------------
// Prepare OCC geometry
// -------------------------
TopoDS_Shape _shape=...;
netgen::OCCGeometry occgeo;
occgeo.shape = _shape;
occgeo.changed = 1;
occgeo.BuildFMap();
// the shape need to have a starting mesh
if (!BRepTools::Triangulation(_shape,1)) {
BRepMesh::Mesh(_shape,1.);
}

Bnd_Box bb;
BRepBndLib::Add (_shape, bb);
double x1,y1,z1,x2,y2,z2;
bb.Get (x1,y1,z1,x2,y2,z2);
netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
occgeo.boundingbox = netgen::Box<3> (p1,p2);

// -------------------------
// Generate the mesh
// -------------------------
netgen::Mesh *ngMesh = NULL;
// we start always with ANALYSE,
// but end depending on _optimize and _isVolume
int startWith = netgen::MESHCONST_ANALYSE;
int endWith = netgen::MESHCONST_OPTSURFACE;
char *optstr=NULL;
int err = 0;
try {
err = netgen::OCCGenerateMesh(occgeo, ngMesh,netgen::mparam, startWith, endWith);
if (!err && netgen::mparam.secondorder > 0) {
netgen::OCCRefinementSurfaces ref (occgeo);
ref.MakeSecondOrder (*ngMesh);
}
}
catch (netgen::NgException exc) {
std::cout <<"Exception in NETGEN: " << exc.What()<GetNP()>0) {

TColgp_Array1OfPnt pts(1,ngMesh->GetNP());
for (int i=1; i<=ngMesh->GetNP(); i++) {
netgen::MeshPoint pt=ngMesh->Point(i);
pts.SetValue(i,gp_Pnt(pt(0),pt(1),pt(2)));
}

Poly_Array1OfTriangle Tri(1,ngMesh->GetNSE());
for (int i=1; i<=ngMesh->GetNSE(); i++) {
// get surface element
netgen::Element2d elem=ngMesh->SurfaceElement(i);
if (elem.GetNP()>=3) {
Standard_Integer N1=elem[0].GetInt();
Standard_Integer N2=elem[1].GetInt();
Standard_Integer N3=elem[2].GetInt();

Tri.SetValue(i,Poly_Triangle(N1,N2,N3));
}
}
aTriangulation=new Poly_Triangulation(pts,Tri);
}
nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
nglib::Ng_Exit();