GModelIO_OCC.cpp

// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.

#include "GmshConfig.h"
#include "GmshMessage.h"
#include "GModelIO_OCC.h"
#include "Context.h"
#include "OCCVertex.h"
#include "OCCEdge.h"
#include "OCCFace.h"
#include "OCCRegion.h"
#include "MElement.h"
#include "MLine.h"
#include "OpenFile.h"
#include "StringUtils.h"
#include "ExtrudeParams.h"

#if defined(HAVE_OCC)

#include <Bnd_Box.hxx>
#include <BRepAlgoAPI_Common.hxx>
#include <BRepAlgoAPI_Cut.hxx>
#include <BRepAlgoAPI_Fuse.hxx>
#include <BRepAlgoAPI_Section.hxx>
#include <BRepBndLib.hxx>
#include <BRepBuilderAPI_Copy.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepBuilderAPI_MakeShell.hxx>
#include <BRepBuilderAPI_MakeSolid.hxx>
#include <BRepBuilderAPI_MakeVertex.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_Sewing.hxx>
#include <BRepBuilderAPI_Transform.hxx>
#include <BRepBuilderAPI_GTransform.hxx>
#include <BRepCheck_Analyzer.hxx>
#include <BRepFilletAPI_MakeFillet.hxx>
#include <BRepGProp.hxx>
#include <BRepLib.hxx>
#include <BRepMesh_FastDiscret.hxx>
#include <BRepOffsetAPI_MakeFilling.hxx>
#include <BRepOffsetAPI_MakePipe.hxx>
#include <BRepOffsetAPI_MakeThickSolid.hxx>
#include <BRepOffsetAPI_Sewing.hxx>
#include <BRepOffsetAPI_ThruSections.hxx>
#include <BRepPrimAPI_MakeBox.hxx>
#include <BRepPrimAPI_MakeCone.hxx>
#include <BRepPrimAPI_MakeCylinder.hxx>
#include <BRepPrimAPI_MakePrism.hxx>
#include <BRepPrimAPI_MakeRevol.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>
#include <BRepPrimAPI_MakeTorus.hxx>
#include <BRepPrimAPI_MakeWedge.hxx>
#include <BRepTools.hxx>
#include <BRep_Tool.hxx>
#include <ElCLib.hxx>
#include <GProp_GProps.hxx>
#include <GeomAPI_PointsToBSpline.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_BezierCurve.hxx>
#include <Geom_Circle.hxx>
#include <Geom_Ellipse.hxx>
#include <Geom_Surface.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <IGESControl_Reader.hxx>
#include <IGESControl_Writer.hxx>
#include <Poly_Triangulation.hxx>
#include <Poly_Triangle.hxx>
#include <STEPControl_Reader.hxx>
#include <STEPControl_Writer.hxx>
#include <ShapeBuild_ReShape.hxx>
#include <ShapeFix_FixSmallFace.hxx>
#include <ShapeFix_Shape.hxx>
#include <ShapeFix_Wireframe.hxx>
#include <Standard_Version.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopTools_DataMapIteratorOfDataMapOfIntegerShape.hxx>
#include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopoDS.hxx>
#include <gce_MakeCirc.hxx>
#include <gce_MakeElips.hxx>
#include <gce_MakePln.hxx>

#if OCC_VERSION_HEX < 0x060900
#error "Gmsh requires OpenCASCADE >= 6.9"
#endif

OCC_Internals::OCC_Internals()
{
  for(int i = 0; i < 6; i++) _maxTag[i] = 0;
  _changed = true;
}

void OCC_Internals::reset()
{
  for(int i = 0; i < 6; i++) _maxTag[i] = 0;
  _meshAttr.Clear();
  _somap.Clear(); _shmap.Clear(); _fmap.Clear(); _wmap.Clear(); _emap.Clear();
  _vmap.Clear();
  _vertexTag.Clear(); _edgeTag.Clear(); _faceTag.Clear(); _solidTag.Clear();
  _tagVertex.Clear(); _tagEdge.Clear(); _tagFace.Clear(); _tagSolid.Clear();
  _wireTag.Clear(); _shellTag.Clear();
  _tagWire.Clear(); _tagShell.Clear();
  _changed = true;
}

void OCC_Internals::setMaxTag(int dim, int val)
{
  if(dim < -2 || dim > 3) return;
  _maxTag[dim + 2] = val;
}

int OCC_Internals::getMaxTag(int dim) const
{
  if(dim < -2 || dim > 3) return 0;
  return _maxTag[dim + 2];
}

void OCC_Internals::_recomputeMaxTag(int dim)
{
  if(dim < -2 || dim > 3) return;
  _maxTag[dim + 2] = 0;
  TopTools_DataMapIteratorOfDataMapOfIntegerShape exp;
  switch(dim){
  case 0: exp.Initialize(_tagVertex); break;
  case 1: exp.Initialize(_tagEdge); break;
  case 2: exp.Initialize(_tagFace); break;
  case 3: exp.Initialize(_tagSolid); break;
  case -1: exp.Initialize(_tagWire); break;
  case -2: exp.Initialize(_tagShell); break;
  }
  for(; exp.More(); exp.Next())
    _maxTag[dim + 2] = std::max(_maxTag[dim + 2], exp.Key());
}

void OCC_Internals::bind(TopoDS_Vertex vertex, int tag, bool recursive)
{
  if(_vertexTag.IsBound(vertex) && _vertexTag.Find(vertex) != tag){
    Msg::Debug("OpenCASCADE vertex %d is already bound to another tag", tag);
  }
  else{
    _vertexTag.Bind(vertex, tag);
    _tagVertex.Bind(tag, vertex);
    setMaxTag(0, tag);
    _changed = true;
  }
}

void OCC_Internals::bind(TopoDS_Edge edge, int tag, bool recursive)
{
  if(_edgeTag.IsBound(edge) && _edgeTag.Find(edge) != tag){
    Msg::Debug("OpenCASCADE edge %d is already bound to another tag", tag);
  }
  else{
    _edgeTag.Bind(edge, tag);
    _tagEdge.Bind(tag, edge);
    setMaxTag(1, tag);
    _changed = true;
  }
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(edge, TopAbs_VERTEX); exp0.More(); exp0.Next()){
      TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
      if(!_vertexTag.IsBound(vertex)){
        int t = getMaxTag(0) + 1;
        bind(vertex, t, recursive);
      }
    }
  }
}

void OCC_Internals::bind(TopoDS_Wire wire, int tag, bool recursive)
{
  if(_wireTag.IsBound(wire) && _wireTag.Find(wire) != tag){
    Msg::Debug("OpenCASCADE wire %d is already bound to anthor tag", tag);
  }
  else{
    _wireTag.Bind(wire, tag);
    _tagWire.Bind(tag, wire);
    setMaxTag(-1, tag);
    _changed = true;
  }
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(wire, TopAbs_EDGE); exp0.More(); exp0.Next()){
      TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
      if(!_edgeTag.IsBound(edge)){
        int t = getMaxTag(1) + 1;
        bind(edge, t, recursive);
      }
    }
  }
}

void OCC_Internals::bind(TopoDS_Face face, int tag, bool recursive)
{
  if(_faceTag.IsBound(face) && _faceTag.Find(face) != tag){
    Msg::Debug("OpenCASCADE face %d is already bound to another tag", tag);
  }
  else{
    _faceTag.Bind(face, tag);
    _tagFace.Bind(tag, face);
    setMaxTag(2, tag);
    _changed = true;
  }
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(face, TopAbs_WIRE); exp0.More(); exp0.Next()){
      TopoDS_Wire wire = TopoDS::Wire(exp0.Current());
      if(!_wireTag.IsBound(wire)){
        int t = getMaxTag(-1) + 1;
        bind(wire, t, recursive);
      }
    }
    for(exp0.Init(face, TopAbs_EDGE); exp0.More(); exp0.Next()){
      TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
      if(!_edgeTag.IsBound(edge)){
        int t = getMaxTag(1) + 1;
        bind(edge, t, recursive);
      }
    }
  }
}

void OCC_Internals::bind(TopoDS_Shell shell, int tag, bool recursive)
{
  if(_shellTag.IsBound(shell) && _shellTag.Find(shell) != tag){
    Msg::Debug("OpenCASCADE shell %d is already bound to another tag", tag);
  }
  else{
    _shellTag.Bind(shell, tag);
    _tagShell.Bind(tag, shell);
    setMaxTag(-2, tag);
    _changed = true;
  }
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(shell, TopAbs_FACE); exp0.More(); exp0.Next()){
      TopoDS_Face face = TopoDS::Face(exp0.Current());
      if(!_faceTag.IsBound(face)){
        int t = getMaxTag(2) + 1;
        bind(face, t, recursive);
      }
    }
  }
}

void OCC_Internals::bind(TopoDS_Solid solid, int tag, bool recursive)
{
  if(_solidTag.IsBound(solid) && _solidTag.Find(solid) != tag){
    Msg::Debug("OpenCASCADE solid %d is already bound to another tag", tag);
  }
  else{
    _solidTag.Bind(solid, tag);
    _tagSolid.Bind(tag, solid);
    setMaxTag(3, tag);
    _changed = true;
  }
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(solid, TopAbs_SHELL); exp0.More(); exp0.Next()){
      TopoDS_Shell shell = TopoDS::Shell(exp0.Current());
      if(!_shellTag.IsBound(shell)){
        int t = getMaxTag(-2) + 1;
        bind(shell, t, recursive);
      }
    }
    for(exp0.Init(solid, TopAbs_FACE); exp0.More(); exp0.Next()){
      TopoDS_Face face = TopoDS::Face(exp0.Current());
      if(!_faceTag.IsBound(face)){
        int t = getMaxTag(3) + 1;
        bind(face, t, recursive);
      }
    }
  }
}

void OCC_Internals::bind(TopoDS_Shape shape, int dim, int tag, bool recursive)
{
  switch(dim){
  case 0: bind(TopoDS::Vertex(shape), tag, recursive); break;
  case 1: bind(TopoDS::Edge(shape), tag, recursive); break;
  case 2: bind(TopoDS::Face(shape), tag, recursive); break;
  case 3: bind(TopoDS::Solid(shape), tag, recursive); break;
  case -1: bind(TopoDS::Wire(shape), tag, recursive); break;
  case -2: bind(TopoDS::Shell(shape), tag, recursive); break;
  default: break;
  }
}

void OCC_Internals::unbind(TopoDS_Vertex vertex, int tag, bool recursive)
{
  TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagEdge);
  for(; exp0.More(); exp0.Next()){
    TopoDS_Edge edge = TopoDS::Edge(exp0.Value());
    TopExp_Explorer exp1;
    for(exp1.Init(edge, TopAbs_VERTEX); exp1.More(); exp1.Next()){
      if(exp1.Current().IsSame(vertex)) return;
    }
  }
  _vertexTag.UnBind(vertex);
  _tagVertex.UnBind(tag);
  _toRemove.insert(std::pair<int, int>(0, tag));
  _recomputeMaxTag(0);
  _changed = true;
}

void OCC_Internals::unbind(TopoDS_Edge edge, int tag, bool recursive)
{
  TopTools_DataMapIteratorOfDataMapOfIntegerShape exp2(_tagFace);
  for(; exp2.More(); exp2.Next()){
    TopoDS_Face face = TopoDS::Face(exp2.Value());
    TopExp_Explorer exp1;
    for(exp1.Init(face, TopAbs_EDGE); exp1.More(); exp1.Next()){
      if(exp1.Current().IsSame(edge)) return;
    }
  }
  _edgeTag.UnBind(edge);
  _tagEdge.UnBind(tag);
  _toRemove.insert(std::pair<int, int>(1, tag));
  _recomputeMaxTag(1);
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(edge, TopAbs_VERTEX); exp0.More(); exp0.Next()){
      TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
      if(_vertexTag.IsBound(vertex)){
        int t = _vertexTag.Find(vertex);
        unbind(vertex, t, recursive);
      }
    }
  }
  _changed = true;
}

void OCC_Internals::unbind(TopoDS_Wire wire, int tag, bool recursive)
{
  TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagFace);
  for(; exp0.More(); exp0.Next()){
    TopoDS_Face face = TopoDS::Face(exp0.Value());
    TopExp_Explorer exp1;
    for(exp1.Init(face, TopAbs_WIRE); exp1.More(); exp1.Next()){
      if(exp1.Current().IsSame(wire)) return;
    }
  }
  _wireTag.UnBind(wire);
  _tagWire.UnBind(tag);
  _toRemove.insert(std::pair<int, int>(-1, tag));
  _recomputeMaxTag(-1);
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(wire, TopAbs_EDGE); exp0.More(); exp0.Next()){
      TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
      if(_edgeTag.IsBound(edge)){
        int t = _edgeTag.Find(edge);
        unbind(edge, t, recursive);
      }
    }
  }
  _changed = true;
}

void OCC_Internals::unbind(TopoDS_Face face, int tag, bool recursive)
{
  TopTools_DataMapIteratorOfDataMapOfIntegerShape exp2(_tagSolid);
  for(; exp2.More(); exp2.Next()){
    TopoDS_Solid solid = TopoDS::Solid(exp2.Value());
    TopExp_Explorer exp1;
    for(exp1.Init(solid, TopAbs_FACE); exp1.More(); exp1.Next()){
      if(exp1.Current().IsSame(face)) return;
    }
  }
  _faceTag.UnBind(face);
  _tagFace.UnBind(tag);
  _toRemove.insert(std::pair<int, int>(2, tag));
  _recomputeMaxTag(2);
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(face, TopAbs_WIRE); exp0.More(); exp0.Next()){
      TopoDS_Wire wire = TopoDS::Wire(exp0.Current());
      if(_wireTag.IsBound(wire)){
        int t = _wireTag.Find(wire);
        unbind(wire, t, recursive);
      }
    }
    for(exp0.Init(face, TopAbs_EDGE); exp0.More(); exp0.Next()){
      TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
      if(_edgeTag.IsBound(edge)){
        int t = _edgeTag.Find(edge);
        unbind(edge, t, recursive);
      }
    }
  }
  _changed = true;
}

void OCC_Internals::unbind(TopoDS_Shell shell, int tag, bool recursive)
{
  TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagSolid);
  for(; exp0.More(); exp0.Next()){
    TopoDS_Solid solid = TopoDS::Solid(exp0.Value());
    TopExp_Explorer exp1;
    for(exp1.Init(solid, TopAbs_SHELL); exp1.More(); exp1.Next()){
      if(exp1.Current().IsSame(shell)) return;
    }
  }
  _shellTag.UnBind(shell);
  _tagShell.UnBind(tag);
  _toRemove.insert(std::pair<int, int>(-2, tag));
  _recomputeMaxTag(-2);
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(shell, TopAbs_FACE); exp0.More(); exp0.Next()){
      TopoDS_Face face = TopoDS::Face(exp0.Current());
      if(_faceTag.IsBound(face)){
        int t = _faceTag.Find(face);
        unbind(face, t, recursive);
      }
    }
  }
  _changed = true;
}

void OCC_Internals::unbind(TopoDS_Solid solid, int tag, bool recursive)
{
  _solidTag.UnBind(solid);
  _tagSolid.UnBind(tag);
  _toRemove.insert(std::pair<int, int>(3, tag));
  _recomputeMaxTag(3);
  if(recursive){
    TopExp_Explorer exp0;
    for(exp0.Init(solid, TopAbs_SHELL); exp0.More(); exp0.Next()){
      TopoDS_Shell shell = TopoDS::Shell(exp0.Current());
      if(_shellTag.IsBound(shell)){
        int t = _shellTag.Find(shell);
        unbind(shell, t, recursive);
      }
    }
    for(exp0.Init(solid, TopAbs_FACE); exp0.More(); exp0.Next()){
      TopoDS_Face face = TopoDS::Face(exp0.Current());
      if(_faceTag.IsBound(face)){
        int t = _faceTag.Find(face);
        unbind(face, t, recursive);
      }
    }
  }
  _changed = true;
}

void OCC_Internals::unbind(TopoDS_Shape shape, int dim, int tag, bool recursive)
{
  switch(dim){
  case 0: unbind(TopoDS::Vertex(shape), tag, recursive); break;
  case 1: unbind(TopoDS::Edge(shape), tag, recursive); break;
  case 2: unbind(TopoDS::Face(shape), tag, recursive); break;
  case 3: unbind(TopoDS::Solid(shape), tag, recursive); break;
  case -1: unbind(TopoDS::Wire(shape), tag, recursive); break;
  case -2: unbind(TopoDS::Shell(shape), tag, recursive); break;
  default: break;
  }
}

void OCC_Internals::_multiBind(TopoDS_Shape shape, int tag,
                               std::vector<std::pair<int, int> > &outDimTags,
                               bool returnHighestDimOnly, bool recursive,
                               bool returnNewOnly)
{
  TopExp_Explorer exp0;
  int count = 0;
  for(exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()){
    TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
    bool exists = false;
    int t = tag;
    if(t <= 0){
      if(_solidTag.IsBound(solid)){
        t = _solidTag.Find(solid);
        exists = true;
      }
      else
        t = getMaxTag(3) + 1;
    }
    else if(count){
      Msg::Error("Cannot bind multiple regions to single tag %d", t);
      return;
    }
    if(!exists)
      bind(solid, t, recursive);
    if(!exists || !returnNewOnly)
      outDimTags.push_back(std::pair<int, int>(3, t));
    count++;
  }
  if(returnHighestDimOnly && count) return;
  for(exp0.Init(shape, TopAbs_FACE); exp0.More(); exp0.Next()){
    TopoDS_Face face = TopoDS::Face(exp0.Current());
    bool exists = false;
    int t = tag;
    if(t <= 0){
      if(_faceTag.IsBound(face)){
        t = _faceTag.Find(face);
        exists = true;
      }
      else
        t = getMaxTag(2) + 1;
    }
    else if(count){
      Msg::Error("Cannot bind multiple faces to single tag %d", t);
      return;
    }
    if(!exists)
      bind(face, t, recursive);
    if(!exists || !returnNewOnly)
      outDimTags.push_back(std::pair<int, int>(2, t));
    count++;
  }
  if(returnHighestDimOnly && count) return;
  for(exp0.Init(shape, TopAbs_EDGE); exp0.More(); exp0.Next()){
    TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
    bool exists = false;
    int t = tag;
    if(t <= 0){
      if(_edgeTag.IsBound(edge)){
        t = _edgeTag.Find(edge);
        exists = true;
      }
      else
        t = getMaxTag(1) + 1;
    }
    else if(count){
      Msg::Error("Cannot bind multiple edges to single tag %d", t);
      return;
    }
    if(!exists)
      bind(edge, t, recursive);
    if(!exists || !returnNewOnly)
      outDimTags.push_back(std::pair<int, int>(1, t));
    count++;
  }
  if(returnHighestDimOnly && count) return;
  for(exp0.Init(shape, TopAbs_VERTEX); exp0.More(); exp0.Next()){
    TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
    bool exists = false;
    int t = tag;
    if(t <= 0){
      if(_vertexTag.IsBound(vertex)){
        t = _vertexTag.Find(vertex);
        exists = false;
      }
      t = getMaxTag(0) + 1;
    }
    else if(count){
      Msg::Error("Cannot bind multiple vertices to single tag %d", t);
      return;
    }
    if(!exists)
      bind(vertex, t, recursive);
    if(!exists || !returnNewOnly)
      outDimTags.push_back(std::pair<int, int>(0, t));
    count++;
  }
}

bool OCC_Internals::_isBound(int dim, int tag)
{
  switch(dim){
  case 0 : return _tagVertex.IsBound(tag);
  case 1 : return _tagEdge.IsBound(tag);
  case 2 : return _tagFace.IsBound(tag);
  case 3 : return _tagSolid.IsBound(tag);
  case -1 : return _tagWire.IsBound(tag);
  case -2 : return _tagShell.IsBound(tag);
  default: return false;
  }
}

bool OCC_Internals::_isBound(int dim, TopoDS_Shape shape)
{
  switch(dim){
  case 0 : return _vertexTag.IsBound(shape);
  case 1 : return _edgeTag.IsBound(shape);
  case 2 : return _faceTag.IsBound(shape);
  case 3 : return _solidTag.IsBound(shape);
  case -1 : return _wireTag.IsBound(shape);
  case -2 : return _shellTag.IsBound(shape);
  default: return false;
  }
}

TopoDS_Shape OCC_Internals::_find(int dim, int tag)
{
  switch(dim){
  case 0: return _tagVertex.Find(tag);
  case 1: return _tagEdge.Find(tag);
  case 2: return _tagFace.Find(tag);
  case 3: return _tagSolid.Find(tag);
  case -1: return _tagWire.Find(tag);
  case -2: return _tagShell.Find(tag);
  default: return TopoDS_Shape();
  }
}

int OCC_Internals::_find(int dim, TopoDS_Shape shape)
{
  switch(dim){
  case 0 : return _vertexTag.Find(shape);
  case 1 : return _edgeTag.Find(shape);
  case 2 : return _faceTag.Find(shape);
  case 3 : return _solidTag.Find(shape);
  case -1 : return _wireTag.Find(shape);
  case -2 : return _shellTag.Find(shape);
  default : return -1;
  }
}

bool OCC_Internals::addVertex(int &tag, double x, double y, double z,
                              double meshSize)
{
  if(tag >= 0 && _tagVertex.IsBound(tag)){
    Msg::Error("OpenCASCADE vertex with tag %d already exists", tag);
    return false;
  }
  TopoDS_Vertex result;
  try{
    gp_Pnt aPnt(x, y, z);
    BRepBuilderAPI_MakeVertex v(aPnt);
    v.Build();
    if(!v.IsDone()){
      Msg::Error("Could not create vertex");
      return false;
    }
    result = v.Vertex();
  }
  catch(Standard_Failure &err){
    Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
    return false;
  }
  if(meshSize > 0 && meshSize < MAX_LC)
    _meshAttr.Bind(result, meshAttr(meshSize));
  if(tag < 0) tag = getMaxTag(0) + 1;
  bind(result, tag, true);
  return true;
}

bool OCC_Internals::addLine(int &tag, int startTag, int endTag)
{
  if(tag >= 0 && _tagEdge.IsBound(tag)){
    Msg::Error("OpenCASCADE edge with tag %d already exists", tag);
    return false;
  }
  if(!_tagVertex.IsBound(startTag)){
    Msg::Error("Unknown OpenCASCADE vertex with tag %d", startTag);
    return false;
  }
  if(!_tagVertex.IsBound(endTag)){
    Msg::Error("Unknown OpenCASCADE vertex with tag %d", endTag);
    return false;
  }

  TopoDS_Edge result;
  try{
    TopoDS_Vertex start = TopoDS::Vertex(_tagVertex.Find(startTag));
    TopoDS_Vertex end = TopoDS::Vertex(_tagVertex.Find(endTag));
    BRepBuilderAPI_MakeEdge e(start, end);
    e.Build();
    if(!e.IsDone()){
      Msg::Error("Could not create edge");
      return false;
    }
    result = e.Edge();
  }
  catch(Standard_Failure &err){
    Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
    return false;
  }
  if(tag < 0) tag = getMaxTag(1) + 1;
  bind(result, tag, true);
  return true;
}

bool OCC_Internals::addLine(int &tag, const std::vector<int> &vertexTags)
{
  if(vertexTags.size() == 2)
    return addLine(tag, vertexTags[0], vertexTags[1]);

  // FIXME: if tag < 0 we could create multiple lines
  Msg::Error("OpenCASCADE polyline currently not supported");
  return false;
}

bool OCC_Internals::_addArc(int &tag, int startTag, int centerTag, int endTag,
                            int mode)
{
  if(tag >= 0 && _tagEdge.IsBound(tag)){
    Msg::Error("OpenCASCADE edge with tag %d already exists", tag);
    return false;
  }
  if(!_tagVertex.IsBound(startTag)){
    Msg::Error("Unknown OpenCASCADE vertex with tag %d", startTag);
    return false;
  }
  if(!_tagVertex.IsBound(centerTag)){
    Msg::Error("Unknown OpenCASCADE vertex with tag %d", centerTag);
    return false;
  }
  if(!_tagVertex.IsBound(endTag)){
    Msg::Error("Unknown OpenCASCADE vertex with tag %d", endTag);
    return false;
  }

  TopoDS_Edge result;
  try{
    TopoDS_Vertex start = TopoDS::Vertex(_tagVertex.Find(startTag));
    TopoDS_Vertex center = TopoDS::Vertex(_tagVertex.Find(centerTag));
    TopoDS_Vertex end = TopoDS::Vertex(_tagVertex.Find(endTag));
    gp_Pnt aP1 = BRep_Tool::Pnt(start);
    gp_Pnt aP2 = BRep_Tool::Pnt(center);
    gp_Pnt aP3 = BRep_Tool::Pnt(end);
    Handle(Geom_TrimmedCurve) arc;
    if(mode == 0){ // circle
      Standard_Real Radius = aP1.Distance(aP2);
      gce_MakeCirc MC(aP2, gce_MakePln(aP1, aP2, aP3).Value(), Radius);
      if(!MC.IsDone()){
        Msg::Error("Could not build circle");
        return false;
      }
      const gp_Circ &Circ = MC.Value();
      Standard_Real Alpha1 = ElCLib::Parameter(Circ, aP1);
      Standard_Real Alpha2 = ElCLib::Parameter(Circ, aP3);
      Handle(Geom_Circle) C = new Geom_Circle(Circ);
      arc = new Geom_TrimmedCurve(C, Alpha1, Alpha2, false);
    }
    else{
      gce_MakeElips ME(aP1, aP3, aP2);
      if(!ME.IsDone()){
        Msg::Error("Could not build ellipse");
        return false;
      }
      const gp_Elips &Elips = ME.Value();
      Standard_Real Alpha1 = ElCLib::Parameter(Elips, aP1);
      Standard_Real Alpha2 = ElCLib::Parameter(Elips, aP3);
      Handle(Geom_Ellipse) E = new Geom_Ellipse(Elips);
      arc = new Geom_TrimmedCurve(E, Alpha1, Alpha2, true);
    }
    BRepBuilderAPI_MakeEdge e(arc, start, end);
    e.Build();
    if(!e.IsDone()){
      Msg::Error("Could not create %s arc", mode ? "ellipse" : "circle");
      return false;
    }
    result = e.Edge();
  }
  catch(Standard_Failure &err){
    Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
    return false;
  }
  if(tag < 0) tag = getMaxTag(1) + 1;
  bind(result, tag, true);
  return true;
}

bool OCC_Internals::addCircleArc(int &tag, int startTag, int centerTag, int endTag)
{
  return _addArc(tag, startTag, centerTag, endTag, 0);
}

bool OCC_Internals::addEllipseArc(int &tag, int startTag, int centerTag, int endTag)
{
  return _addArc(tag, startTag, centerTag, endTag, 1);
}

bool OCC_Internals::addCircle(int &tag, double x, double y, double z, double r,
                              double angle1, double angle2)
{
  if(tag >= 0 && _tagEdge.IsBound(tag)){
    Msg::Error("OpenCASCADE edge with tag %d already exists", tag);
    return false;
  }

  TopoDS_Edge result;
  try{
    gp_Dir N_dir(0., 0., 1.);
    gp_Dir x_dir(1., 0., 0.);
    gp_Pnt center(x, y, z);
    gp_Ax2 axis(center, N_dir, x_dir);
    gp_Circ circ(axis, r);
    if(angle1 == 0. && angle2 == 2 * M_PI){
      result = BRepBuilderAPI_MakeEdge(circ);
    }
    else{
      Handle(Geom_Circle) C = new Geom_Circle(circ);
      Handle(Geom_TrimmedCurve) arc = new Geom_TrimmedCurve(C, angle1, angle2, false);
      BRepBuilderAPI_MakeEdge e(arc);
      if(!e.IsDone()){
        Msg::Error("Could not create circle arc");
        return false;
      }
      result = e.Edge();
    }
  }
  catch(Standard_Failure &err){
    Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
    return false;
  }
  if(tag < 0) tag = getMaxTag(1) + 1;
  bind(result, tag, true);
  return true;
}

bool OCC_Internals::addEllipse(int &tag, double x, double y, double z, double rx,
                               double ry, double angle1, double angle2)
{
  if(tag >= 0 && _tagEdge.IsBound(tag)){
    Msg::Error("OpenCASCADE edge with tag %d already exists", tag);
    return false;
  }
  if(ry > rx){
    Msg::Error("Major radius rx should be larger than minor radius ry");
    return false;
  }

  TopoDS_Edge result;
  try{
    gp_Dir N_dir(0., 0., 1.);
    gp_Dir x_dir(1., 0., 0.);
    gp_Pnt center(x, y, z);
    gp_Ax2 axis(center, N_dir, x_dir);
    gp_Elips elips(axis, rx, ry);
    if(angle1 == 0 && angle2 == 2 * M_PI){
      result = BRepBuilderAPI_MakeEdge(elips);
    }
    else{
      Handle(Geom_Ellipse) E = new Geom_Ellipse(elips);
      Handle(Geom_TrimmedCurve) arc = new Geom_TrimmedCurve(E, angle1, angle2, true);
      BRepBuilderAPI_MakeEdge e(arc);
      if(!e.IsDone()){
        Msg::Error("Could not create ellipse arc");
        return false;
      }
      result = e.Edge();
    }
  }
  catch(Standard_Failure &err){
    Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
    return false;
  }
  if(tag < 0) tag = getMaxTag(1) + 1;
  bind(result, tag, true);
  return true;
}

bool OCC_Internals::_addSpline(int &tag, const std::vector<int> &vertexTags, int mode)
{
  if(tag >= 0 && _tagEdge.IsBound(tag)){
    Msg::Error("OpenCASCADE edge with tag %d already exists", tag);
    return false;
  }
  if(vertexTags.size() < 2 || vertexTags.size() > 20){
    Msg::Error("Number of control points should be in [2,20]");
    return false;
  }

  TopoDS_Edge result;
  try{
    TColgp_Array1OfPnt ctrlPoints(1, vertexTags.size());
    TopoDS_Vertex start, end;
    for(unsigned int i = 0; i < vertexTags.size(); i++){
      if(!_tagVertex.IsBound(vertexTags[i])){
        Msg::Error("Unknown OpenCASCADE vertex with tag %d", vertexTags[i]);
        return false;
      }
      TopoDS_Vertex vertex = TopoDS::Vertex(_tagVertex.Find(vertexTags[i]));
      ctrlPoints.SetValue(i + 1, BRep_Tool::Pnt(vertex));
      if(i == 0) start = vertex;
      if(i == vertexTags.size() - 1) end = vertex;
    }
    if(mode == 0){ // BSpline through points (called "Spline" in Gmsh)
      Handle(Geom_BSplineCurve) curve = GeomAPI_PointsToBSpline(ctrlPoints).Curve();
      BRepBuilderAPI_MakeEdge e(curve, start, end);
      if(!e.IsDone()){
        Msg::Error("Could not create spline");
        return false;
      }
      result = e.Edge();
    }
    else if(mode == 1){
      Handle(Geom_BezierCurve) curve = new Geom_BezierCurve(ctrlPoints);
      BRepBuilderAPI_MakeEdge e(curve, start, end);
      if(!e.IsDone()){
        Msg::Error("Could not create Bezier curve");
        return false;
      }
      result = e.Edge();
    }
    else if(mode == 2){
      // TODO: BSpline treat periodic case, allow to set order, etc.
      // Handle(Geom_BSplineCurve) curve = new Geom_BSplineCurve(ctrlPoints, ...);
      // ...
      Msg::Error("OpenCASCADE BSpline not implemented yet");
      return false;
    }
  }
  catch(Standard_Failure &err){
    Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
    return false;
  }
  if(tag < 0) tag = getMaxTag(1) + 1;
  bind(result, tag, true);
  return true;
}

bool OCC_Internals::addSpline(int &tag, const std::vector<int> &vertexTags)
{
  return _addSpline(tag, vertexTags, 0);
}

bool OCC_Internals::addBezier(int &tag, const std::vector<int> &vertexTags)
{
  return _addSpline(tag, vertexTags, 1);
}

bool OCC_Internals::addBSpline(int &tag, const std::vector<int> &vertexTags)
{
  return _addSpline(tag, vertexTags, 2);
}

bool OCC_Internals::addWire(int &tag, const std::vector<int> &edgeTags,
                            bool checkClosed)
{
  if(tag >= 0 && _tagWire.IsBound(tag)){
    Msg::Error("OpenCASCADE wire or line loop with tag %d already exists", tag);
    return false;
  }

  TopoDS_Wire result;
  try{
    BRepBuilderAPI_MakeWire w;
    for (unsigned int i = 0; i < edgeTags.size(); i++) {
      if(!_tagEdge.IsBound(edgeTags[i])){
        Msg::Error("Unknown OpenCASCADE edge with tag %d", edgeTags[i]);
        return false;
      }
      TopoDS_Edge edge = TopoDS::Edge(_tagEdge.Find(edgeTags[i]));
      w.Add(edge);
    }
    result = w.Wire();
    if(checkClosed && !result.Closed()){
      Msg::Error("Line Loop is not closed");
      return false;
    }
  }
  catch(Standard_Failure &err){
    Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
    return false;
  }
  if(tag < 0) tag = getMaxTag(-1) + 1;
  bind(result, tag, true);
  return true;
}

bool OCC_Internals::addLineLoop(int &tag, const std::vector<int> &edgeTags)
{
  return addWire(tag, edgeTags, true);
}

bool OCC_Internals::_makeRectangle(TopoDS_Face &result, double x, double y, double z,
                                   double dx, double dy, double roundedRadius)
{
  try{
    TopoDS_Wire wire;
    double x1 = x, y1 = y, z1 = z, x2 = x1 + dx, y2 = y1 + dy;
    if(roundedRadius <= 0.){
      TopoDS_Vertex v1 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y1, z1));
      TopoDS_Vertex v2 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y1, z1));
      TopoDS_Vertex v3 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y2, z1));
      TopoDS_Vertex v4 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y2, z1));
      TopoDS_Edge e1 = BRepBuilderAPI_MakeEdge(v1, v2);
      TopoDS_Edge e2 = BRepBuilderAPI_MakeEdge(v2, v3);
      TopoDS_Edge e3 = BRepBuilderAPI_MakeEdge(v3, v4);
      TopoDS_Edge e4 = BRepBuilderAPI_MakeEdge(v4, v1);
      wire = BRepBuilderAPI_MakeWire(e1, e2, e3, e4);
    }
    else{
      double r = roundedRadius;
      TopoDS_Vertex v1 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1 + r, y1, z1));
      TopoDS_Vertex v2 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2 - r, y1, z1));
      TopoDS_Vertex v3 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y1 + r, z1));
      TopoDS_Vertex v4 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y2 - r, z1));
      TopoDS_Vertex v5 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2 - r, y2, z1));
      TopoDS_Vertex v6 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1 + r, y2, z1));
      TopoDS_Vertex v7 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y2 - r, z1));
      TopoDS_Vertex v8 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y1 + r, z1));
      TopoDS_Edge e1 = BRepBuilderAPI_MakeEdge(v1, v2);
      TopoDS_Edge e2 = BRepBuilderAPI_MakeEdge(v3, v4);
      TopoDS_Edge e3 = BRepBuilderAPI_MakeEdge(v5, v6);
      TopoDS_Edge e4 = BRepBuilderAPI_MakeEdge(v7, v8);
      gp_Pnt c1(x1 + r, y1 + r, z1);