From ca9261dc4d4f3f053b56643b855920148a4c7435 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Sat, 1 May 2010 15:28:11 +0000
Subject: [PATCH] new modeler taks form ...

---
 CMakeLists.txt             |  15 ++
 Common/GmshConfig.h.in     |   1 +
 Common/OpenFile.cpp        |   3 +
 Geo/ACISEdge.cpp           | 107 ++++++++++++
 Geo/ACISEdge.h             |  48 ++++++
 Geo/ACISVertex.cpp         |  56 +++++++
 Geo/ACISVertex.h           |  38 +++++
 Geo/CMakeLists.txt         |   3 +
 Geo/GEntity.h              |   3 +-
 Geo/GModel.cpp             |  28 +++-
 Geo/GModel.h               |  23 ++-
 Geo/GModelFactory.cpp      | 323 +++++++++++++++++++++++++++----------
 Geo/GModelFactory.h        | 104 ++++++------
 Geo/GModelIO_ACIS.cpp      | 149 +++++++++++++++++
 Geo/GModelIO_ACIS.h        |   0
 Geo/GModelIO_OCC.cpp       | 168 +++++--------------
 Geo/GModelIO_OCC.h         |  22 ++-
 Geo/GVertex.cpp            |   7 +
 Geo/OCCIncludes.h          |   1 -
 Geo/OCCRegion.cpp          |  13 ++
 Geo/OCCRegion.h            |   1 +
 Solver/linearSystemPETSc.h |   4 +-
 22 files changed, 843 insertions(+), 274 deletions(-)
 create mode 100644 Geo/ACISEdge.cpp
 create mode 100644 Geo/ACISEdge.h
 create mode 100644 Geo/ACISVertex.cpp
 create mode 100644 Geo/ACISVertex.h
 create mode 100644 Geo/GModelIO_ACIS.cpp
 create mode 100644 Geo/GModelIO_ACIS.h

diff --git a/CMakeLists.txt b/CMakeLists.txt
index ca602e655a..09c0e84771 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -39,6 +39,7 @@ option(ENABLE_MSVC_STATIC_RUNTIME "Use static Visual C++ runtime" OFF)
 option(ENABLE_NATIVE_FILE_CHOOSER "Enable native file chooser in GUI" ON)
 option(ENABLE_NETGEN "Enable Netgen mesh generator" ON)
 option(ENABLE_OCC "Enable Open CASCADE geometrical models" ON)
+option(ENABLE_ACIS "Enable ACIS geometrical models" ON)
 option(ENABLE_OSMESA "Use OSMesa for offscreen rendering" OFF)
 option(ENABLE_PARSER "Build the GEO file parser" ON)
 option(ENABLE_PETSC "Enable PETSc linear algebra solvers" ON)
@@ -694,6 +695,20 @@ if(ENABLE_OCC)
   endif(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
 endif(ENABLE_OCC)
 
+if(ENABLE_ACIS)
+  find_library(ACIS_LIB SpaACIS PATH_SUFFIXES bin/maci386)
+  if(ACIS_LIB)
+    find_path(ACIS_INC "kernapi.hxx" PATH_SUFFIXES include)
+    if(ACIS_INC)
+      set(HAVE_ACIS TRUE)
+      list(APPEND CONFIG_OPTIONS "Acis")
+      list(APPEND EXTERNAL_LIBRARIES ${ACIS_LIB})
+      list(APPEND EXTERNAL_INCLUDES ${ACIS_INC})
+    endif(ACIS_INC)
+  endif(ACIS_LIB)
+endif(ENABLE_ACIS)
+
+
 if(ENABLE_OSMESA)
   find_library(OSMESA_LIB OSMesa)
   if(OSMESA_LIB)
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index aa1f2f034c..d03dbb4ef9 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -7,6 +7,7 @@
 #define _GMSH_CONFIG_H_
 
 #cmakedefine HAVE_64BIT_SIZE_T
+#cmakedefine HAVE_ACIS
 #cmakedefine HAVE_ANN
 #cmakedefine HAVE_BLAS
 #cmakedefine HAVE_CHACO
diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index cca9362c8f..f51b50ad1c 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -284,6 +284,9 @@ int MergeFile(std::string fileName, bool warnIfMissing)
   else if(ext == ".step" || ext == ".STEP" || ext == ".stp" || ext == ".STP"){
     status = GModel::current()->readOCCSTEP(fileName);
   }
+  else if(ext == ".sat" || ext == ".SAT"){
+    status = GModel::current()->readACISSAT(fileName);
+  }
   else if(ext == ".unv" || ext == ".UNV"){
     status = GModel::current()->readUNV(fileName);
   }
diff --git a/Geo/ACISEdge.cpp b/Geo/ACISEdge.cpp
new file mode 100644
index 0000000000..9104f1e91e
--- /dev/null
+++ b/Geo/ACISEdge.cpp
@@ -0,0 +1,107 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include <limits>
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "GModel.h"
+#include "ACISEdge.h"
+#include "Context.h"
+
+#if defined(HAVE_ACIS)
+
+#include <interval.hxx>
+#include <curve.hxx>
+#include <curdef.hxx>
+
+GEdge *getACISEdgeByNativePtr(GModel *model, EDGE * toFind)
+{
+  GModel::eiter it =model->firstEdge();
+  for (; it !=model->lastEdge(); it++){
+    ACISEdge *ed = dynamic_cast<ACISEdge*>(*it);
+    if (ed){
+      if (toFind == ed->getEDGE()){
+	return *it;
+      }
+    }
+  }
+  return 0;
+}
+
+
+ACISEdge::ACISEdge(GModel *model, EDGE* edge, int num, GVertex *v1, GVertex *v2)
+  : GEdge(model, num, v1, v2), _e(edge)
+{
+  
+  SPAinterval interval = _e->get_param_range (); 
+  s0 = interval.start_pt();
+  s1 = interval.end_pt();
+}
+
+Range<double> ACISEdge::parBounds(int i) const
+{ 
+  return Range<double>(s0, s1);
+}
+
+SPoint2 ACISEdge::reparamOnFace(const GFace *face, double epar, int dir) const
+{
+  throw;
+}
+
+GPoint ACISEdge::closestPoint(const SPoint3 &qp, double &param) const{
+  throw;
+}
+
+
+// True if the edge is a seam for the given face
+bool ACISEdge::isSeam(const GFace *face) const
+{
+  throw;
+}
+
+GPoint ACISEdge::point(double par) const
+{
+  CURVE *c = _e->geometry();
+  const curve &equ = c->equation();
+  if (&equ){
+    SPAposition pos;
+    equ.eval(par,pos);
+    return GPoint(pos.coordinate(0),pos.coordinate(1),pos.coordinate(2),this,par);
+  }
+  Msg::Error("Unable to evaluate a point on a curve");
+  return GPoint(0,0,0,this,par);  
+}
+
+SVector3 ACISEdge::firstDer(double par) const
+{  
+  return SVector3(0,0,0);
+}
+
+GEntity::GeomType ACISEdge::geomType() const
+{
+  return Unknown;
+}
+
+int ACISEdge::minimumMeshSegments() const
+{
+  return GEdge::minimumMeshSegments();
+}
+
+int ACISEdge::minimumDrawSegments() const
+{
+  return GEdge::minimumDrawSegments();
+}
+
+double ACISEdge::curvature(double par) const 
+{
+  const double Crv = 1.e-15;
+  return Crv;
+}
+
+void ACISEdge::writeGEO(FILE *fp)
+{
+}
+
+#endif
diff --git a/Geo/ACISEdge.h b/Geo/ACISEdge.h
new file mode 100644
index 0000000000..ff36fc7808
--- /dev/null
+++ b/Geo/ACISEdge.h
@@ -0,0 +1,48 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _ACIS_EDGE_H_
+#define _ACIS_EDGE_H_
+
+#include "GmshConfig.h"
+#include "GEdge.h"
+#include "GModel.h"
+#include "ACISVertex.h"
+#include "Range.h"
+
+class ACISFace;
+
+#if defined(HAVE_ACIS)
+
+#include <position.hxx>
+#include <edge.hxx>
+
+class ACISEdge : public GEdge {
+  EDGE *_e;
+  double s0, s1;
+ public:
+  ACISEdge(GModel *model, EDGE *e, int num, GVertex *v1, GVertex *v2);
+  virtual ~ACISEdge() {}
+  virtual Range<double> parBounds(int i) const;
+  virtual GeomType geomType() const;
+  virtual bool degenerate(int) const { return false; }
+  virtual GPoint point(double p) const;
+  virtual SVector3 firstDer(double par) const;
+  virtual double curvature (double par) const;
+  virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
+  virtual GPoint closestPoint(const SPoint3 &queryPoint, double &param) const;
+  ModelType getNativeType() const { return AcisModel; }
+  void * getNativePtr() const { return (void*)_e; }
+  virtual int minimumMeshSegments () const;
+  virtual int minimumDrawSegments () const;
+  bool isSeam(const GFace *) const;
+  virtual void writeGEO(FILE *fp);
+  EDGE* getEDGE() const {return _e;}
+};
+GEdge *getACISEdgeByNativePtr(GModel *model, EDGE *toFind);
+
+#endif
+
+#endif
diff --git a/Geo/ACISVertex.cpp b/Geo/ACISVertex.cpp
new file mode 100644
index 0000000000..1c77922e3b
--- /dev/null
+++ b/Geo/ACISVertex.cpp
@@ -0,0 +1,56 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "ACISVertex.h"
+#include "MPoint.h"
+
+#if defined(HAVE_ACIS)
+#include <point.hxx>
+
+ACISVertex::ACISVertex(GModel *m, int num, VERTEX *v) 
+  : GVertex(m, num), _v(v)
+{
+  APOINT *p = _v->geometry();
+  const SPAposition &pos = p->coords();
+  _x = pos.coordinate(0);
+  _y = pos.coordinate(1);
+  _z = pos.coordinate(2);
+
+  mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
+  points.push_back(new MPoint(mesh_vertices.back()));
+}
+
+void ACISVertex::setPosition(GPoint &p)
+{
+  _x = p.x();
+  _y = p.y();
+  _z = p.z();
+  if(mesh_vertices.size()){
+    mesh_vertices[0]->x() = p.x();
+    mesh_vertices[0]->y() = p.y();
+    mesh_vertices[0]->z() = p.z();
+  }
+}
+
+SPoint2 ACISVertex::reparamOnFace(const GFace *gf, int dir) const
+{
+  Msg::Fatal("NOT DONE !!!!!!!!!!!\n");
+}
+
+GVertex *getACISVertexByNativePtr(GModel *model, VERTEX* toFind)
+{
+  GModel::viter it =model->firstVertex();
+  for (; it != model->lastVertex(); it++){
+    ACISVertex *av = dynamic_cast<ACISVertex*>(*it);
+    if (av){
+      if (toFind == av->getVERTEX()){
+	return *it;
+      }
+    }
+  }
+  return 0;
+}
+
+#endif
diff --git a/Geo/ACISVertex.h b/Geo/ACISVertex.h
new file mode 100644
index 0000000000..d242883710
--- /dev/null
+++ b/Geo/ACISVertex.h
@@ -0,0 +1,38 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _ACIS_VERTEX_H_
+#define _ACIS_VERTEX_H_
+
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "GVertex.h"
+
+#if defined(HAVE_ACIS)
+
+#include <vertex.hxx>
+
+class ACISVertex : public GVertex {
+ protected:
+  VERTEX* _v;
+  double _x, _y, _z;
+ public:
+  ACISVertex(GModel *m, int num, VERTEX *_v);
+  virtual ~ACISVertex() {}
+  virtual GPoint point() const { return GPoint(x(), y(), z()); }
+  virtual double x() const { return _x; }
+  virtual double y() const { return _y; }
+  virtual double z() const { return _z; }
+  virtual void setPosition(GPoint &p);
+  ModelType getNativeType() const { return AcisModel; }
+  void * getNativePtr() const { return (void*)_v; }
+  virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
+  VERTEX* getVERTEX() { return _v; }
+};
+GVertex *getACISVertexByNativePtr(GModel *model, VERTEX*);
+
+#endif
+
+#endif
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 77208d2648..e7e4ff7b65 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -11,11 +11,14 @@ set(SRC
     OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp
     discreteEdge.cpp discreteFace.cpp discreteRegion.cpp
     fourierEdge.cpp fourierFace.cpp fourierProjectionFace.cpp
+  ACISVertex.cpp
+  ACISEdge.cpp
   GModel.cpp
   GModelFactory.cpp
     GModelIO_Geo.cpp
     GModelIO_Mesh.cpp
     GModelIO_OCC.cpp
+    GModelIO_ACIS.cpp
       OCC_Connect.cpp
     GModelIO_Fourier.cpp
     GModelIO_CGNS.cpp
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index e717303573..3d5c8ac388 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -65,7 +65,8 @@ class GEntity {
     UnknownModel,
     GmshModel,
     FourierModel,
-    OpenCascadeModel
+    OpenCascadeModel,
+    AcisModel
   };
 
   // all known entity types
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e7b7f1590a..58fcc86ac5 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -42,8 +42,7 @@ int GModel::_current = -1;
 
 GModel::GModel(std::string name)
   : _name(name), _visible(1), _octree(0),
-    _geo_internals(0), _occ_internals(0), _fm_internals(0), _factory(0),
-    _fields(0), _currentMeshEntity(0), normals(0)
+    _geo_internals(0), _occ_internals(0), _acis_internals(0), _fm_internals(0),_factory(0), _fields(0), _currentMeshEntity(0), normals(0)
 {
   partitionSize[0] = 0; partitionSize[1] = 0;
   list.push_back(this);
@@ -1453,11 +1452,21 @@ GEdge *GModel::addCircleArc3Points(double x, double y, double z, GVertex *start,
 }
 
 GEdge *GModel::addBezier(GVertex *start, GVertex *end, 
-                         fullMatrix<double> *controlPoints)
+			 std::vector<std::vector<double> > points) 
 {
   if(_factory) 
     return _factory->addSpline(this, GModelFactory::BEZIER, start, end, 
-                               controlPoints);
+                               points);
+  return 0;
+}
+
+GEdge *GModel::addNURBS(GVertex *start, GVertex *end,
+			std::vector<std::vector<double> > points, 
+			std::vector<double> knots,
+			std::vector<double> weights, 
+			std::vector<int> mult){  
+  if(_factory)
+    return _factory->addNURBS(this, start,end,points,knots,weights, mult);
   return 0;
 }
 
@@ -1738,7 +1747,7 @@ static void glueFacesInRegions(GModel *model,
   }
 }
 
-void GModel::glue(const double &eps)
+void GModel::glue(double eps)
 {
   {
     std::multimap<GVertex*,GVertex*> Unique2Duplicates;
@@ -1821,6 +1830,9 @@ void GModel::registerBindings(binding *b)
   cm->setDescription("create a spline going from v1 to v2 and with some control "
                      "points listed in a fullMatrix(N,3)");
   cm->setArgNames("v1", "v2", "controlPoints", NULL);
+  cm = cb->addMethod("addNURBS", &GModel::addNURBS);
+  cm->setDescription("creates a NURBS curve from v1 to v2 with control Points, knots, weights and multiplicities");
+  cm->setArgNames("v1", "v2", "{{poles}}","{knots}","{weights}","{mult}",NULL);
   cm = cb->addMethod("addCircleArcCenter", &GModel::addCircleArcCenter);
   cm->setDescription("create a circle arc going from v1 to v2 with its center "
                      "at (x,y,z)");
@@ -1865,9 +1877,9 @@ void GModel::registerBindings(binding *b)
                      "one (tool). The third parameter tells if a new model has to "
                      "be created");
   cm->setArgNames("tool","createNewGModel",NULL);
-  cm = cb->addMethod("glue", &GModel::glue);
-  cm->setDescription("glue the geometric model using geometric tolerance eps");
-  cm->setArgNames("eps",NULL);
+  //  cm = cb->addMethod("glue", &GModel::glue);
+  //  cm->setDescription("glue the geometric model using geometric tolerance eps");
+  //  cm->setArgNames("eps",NULL);
   cm = cb->addMethod("setAsCurrent", &GModel::setAsCurrent);
   cm->setDescription("set the model as the current (active) one");
   cm = cb->setConstructor<GModel>();
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 65d2af4850..f97b9854bf 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -23,6 +23,7 @@ class Octree;
 class FM_Internals;
 class GEO_Internals;
 class OCC_Internals;
+class ACIS_Internals;
 class smooth_normals;
 class FieldManager;
 class CGNSOptions;
@@ -68,6 +69,10 @@ class GModel
   OCC_Internals *_occ_internals;
   void _deleteOCCInternals();
 
+  // OpenCascade model internal data
+  ACIS_Internals *_acis_internals;
+  void _deleteACISInternals();
+
   // Fourier model internal data
   FM_Internals *_fm_internals;
   void _createFMInternals();
@@ -149,6 +154,7 @@ class GModel
   GEO_Internals *getGEOInternals(){ return _geo_internals; }
   OCC_Internals *getOCCInternals(){ return _occ_internals; }
   FM_Internals *getFMInternals() { return _fm_internals; }
+  ACIS_Internals *getACISInternals(){ return _acis_internals; }
 
   // access characteristic length (mesh size) fields
   FieldManager *getFields(){ return _fields; }
@@ -346,6 +352,13 @@ class GModel
   // mesh the model
   int mesh(int dimension);
 
+  // glue entities in the model
+  // assume a tolerance eps and merge vertices that are too close, 
+  // then merge edges, faces and regions.
+  // the gluer changes the geometric model, so that some pointers
+  // could become invalid !! I think that using references to some 
+  // tables of pointers for bindings e.g. could be better. FIXME !!
+  void glue (double eps);
   // change the entity creation factory
   void setFactory(std::string name);
 
@@ -354,7 +367,12 @@ class GModel
   GEdge *addLine(GVertex *v1, GVertex *v2);
   GEdge *addCircleArcCenter(double x, double y, double z, GVertex *start, GVertex *end);
   GEdge *addCircleArc3Points(double x, double y, double z, GVertex *start, GVertex *end);
-  GEdge *addBezier(GVertex *start, GVertex *end, fullMatrix<double> *controlPoints);
+  GEdge *addBezier(GVertex *start, GVertex *end, std::vector<std::vector<double> > points);
+  GEdge *addNURBS(GVertex *start, GVertex *end,
+		  std::vector<std::vector<double> > points, 
+		  std::vector<double> knots,
+		  std::vector<double> weights, 
+		  std::vector<int> mult);
   GEntity *revolve(GEntity *e, std::vector<double> p1, std::vector<double> p2, double angle);
   GEntity *extrude(GEntity *e, std::vector<double> p1, std::vector<double> p2);
 
@@ -422,6 +440,9 @@ class GModel
   int writeOCCSTEP(const std::string &name);
   int writeOCCBREP(const std::string &name);
   int importOCCShape(const void *shape);
+  
+  // ACIS Model
+  int readACISSAT(const std::string &name);
 
   // Gmsh mesh file format
   int readMSH(const std::string &name);
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index cc89b6f58e..a5a801d0d3 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -25,6 +25,11 @@
 #include <Geom_Circle.hxx>
 #include <Geom_TrimmedCurve.hxx>
 #include <Geom_BezierCurve.hxx>
+#include <Geom_BSplineCurve.hxx>
+#include <BRepBuilderAPI_Transform.hxx>
+#include <TColgp_HArray1OfPnt.hxx>
+#include <TColStd_HArray1OfReal.hxx>
+#include <TColStd_HArray1OfInteger.hxx>
 
 GVertex *OCCFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
 {
@@ -56,24 +61,30 @@ GEdge *OCCFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
     gp_Pnt p2(end->x(),end->y(),end->z());
     TopoDS_Edge occEdge = BRepBuilderAPI_MakeEdge(p1, p2).Edge();
   }  
-  return gm->_occ_internals->addEdgeToModel(gm, occEdge, start, end);
+  return gm->_occ_internals->addEdgeToModel(gm,occEdge);
 }
 
-GEdge *OCCFactory::addCircleArc(GModel *gm, const arcCreationMethod &method,
-                                GVertex *start, GVertex *end, 
-                                const SPoint3 &aPoint)
-{
+GEdge *OCCFactory::addCircleArc (GModel *gm, const arcCreationMethod &method,
+				 GVertex *start, 
+				 GVertex *end, 
+				 const SPoint3 &aPoint) {
   if (!gm->_occ_internals)
     gm->_occ_internals = new OCC_Internals;
-
+  
   gp_Pnt aP1 (start->x(), start->y(), start->z());
   gp_Pnt aP2 (aPoint.x(), aPoint.y(), aPoint.z());
   gp_Pnt aP3 (end->x(), end->y(), end->z());
+  TopoDS_Edge occEdge;
+
+  OCCVertex *occv1 = dynamic_cast<OCCVertex*>(start);
+  OCCVertex *occv2 = dynamic_cast<OCCVertex*>(end);
 
   if (method == GModelFactory::THREE_POINTS){
     GC_MakeArcOfCircle arc(aP1, aP2, aP3);
-    TopoDS_Edge occEdge = BRepBuilderAPI_MakeEdge(arc).Edge();
-    return gm->_occ_internals->addEdgeToModel(gm, occEdge, start, end);
+    if (occv1 && occv2)
+      occEdge = BRepBuilderAPI_MakeEdge(arc,occv1->getShape(),occv2->getShape()).Edge();    
+    else 
+      occEdge = BRepBuilderAPI_MakeEdge(arc).Edge();
   }
   else if (method == GModelFactory::CENTER_START_END){
     Standard_Real Radius = aP1.Distance(aP2);
@@ -83,41 +94,110 @@ GEdge *OCCFactory::addCircleArc(GModel *gm, const arcCreationMethod &method,
     Standard_Real Alpha2 = ElCLib::Parameter(Circ, aP3);
     Handle(Geom_Circle) C = new Geom_Circle(Circ);
     Handle(Geom_TrimmedCurve) arc = new Geom_TrimmedCurve(C, Alpha1, Alpha2, false);
-    TopoDS_Edge occEdge = BRepBuilderAPI_MakeEdge(arc).Edge();
-    return gm->_occ_internals->addEdgeToModel(gm, occEdge, start, end);
+    if (occv1 && occv2)
+      occEdge = BRepBuilderAPI_MakeEdge(arc,occv1->getShape(),occv2->getShape()).Edge();    
+    else 
+      occEdge = BRepBuilderAPI_MakeEdge(arc).Edge();
   }
-  return 0;
+  return gm->_occ_internals->addEdgeToModel(gm,occEdge);
 }
 
 GEdge *OCCFactory::addSpline(GModel *gm, const splineType &type,
                              GVertex *start, GVertex *end, 
-                             fullMatrix<double> *points)
+			     std::vector<std::vector<double> > points)
 {
   if (!gm->_occ_internals)
     gm->_occ_internals = new OCC_Internals;
 
+  TopoDS_Edge occEdge;
+
+  OCCVertex *occv1 = dynamic_cast<OCCVertex*>(start);
+  OCCVertex *occv2 = dynamic_cast<OCCVertex*>(end);
+
   OCCEdge *occEd = 0;
-  int nbControlPoints = points->size1();
+  int nbControlPoints = points.size();
   TColgp_Array1OfPnt ctrlPoints(1, nbControlPoints + 2);
   int index = 1;
   ctrlPoints.SetValue(index++, gp_Pnt(start->x(), start->y(), start->z()));  
   for (int i = 0; i < nbControlPoints; i++) {
-    gp_Pnt aP((*points)(i, 0), (*points)(i, 1), (*points)(i, 2));
+    gp_Pnt aP(points[i][0],points[i][1],points[i][2]);
     ctrlPoints.SetValue(index++, aP);
   }
   ctrlPoints.SetValue(index++, gp_Pnt(end->x(), end->y(), end->z()));  
   if (type == BEZIER) {
     Handle(Geom_BezierCurve) Bez = new Geom_BezierCurve(ctrlPoints);
-    TopoDS_Edge bez = BRepBuilderAPI_MakeEdge(Bez).Edge();
-    return gm->_occ_internals->addEdgeToModel(gm, bez, start, end);
+    if (occv1 && occv2)
+      occEdge = BRepBuilderAPI_MakeEdge(Bez,occv1->getShape(),occv2->getShape()).Edge();    
+    else
+      occEdge = BRepBuilderAPI_MakeEdge(Bez).Edge();
   } 
-  Msg::Error("Non-bezier splines not implemented yet");
+  return gm->_occ_internals->addEdgeToModel(gm,occEdge);
+}
+
+
+GEdge *OCCFactory::addNURBS(GModel *gm,
+			    GVertex *start, GVertex *end,
+			    std::vector<std::vector<double> > points, 
+			    std::vector<double> knots,
+			    std::vector<double> weights, 
+			    std::vector<int> mult){
+  try{ 
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
+  
+  
+  
+  OCCVertex *occv1 = dynamic_cast<OCCVertex*>(start);
+  OCCVertex *occv2 = dynamic_cast<OCCVertex*>(end);
+  
+  int nbControlPoints = points.size() + 2;
+  TColgp_Array1OfPnt  ctrlPoints(1, nbControlPoints);
+  
+
+  TColStd_Array1OfReal _knots    (1, knots.size());
+  TColStd_Array1OfReal _weights  (1, weights.size());
+  TColStd_Array1OfInteger  _mult     (1, mult.size());
+  
+
+  for (int i = 0; i < knots.size(); i++) {
+    _knots.SetValue(i+1, knots[i]);
+  }
+  for (int i = 0; i < weights.size(); i++) {
+    _weights.SetValue(i+1, weights[i]);
+  }
+  int totKnots = 0;
+  for (int i = 0; i < mult.size(); i++) {
+    _mult.SetValue(i+1, mult[i]);   
+    totKnots += mult[i];
+  }
+
+  const int degree = totKnots - nbControlPoints - 1;
+  printf("creation of a nurbs of degree %d with %d control points\n",degree,nbControlPoints);
+  
+  int index = 1;
+  ctrlPoints.SetValue(index++, gp_Pnt(start->x(), start->y(), start->z()));  
+  for (int i = 0; i < points.size(); i++) {
+    gp_Pnt aP(points[i][0],points[i][1],points[i][2]);
+    ctrlPoints.SetValue(index++, aP);
+  }
+  ctrlPoints.SetValue(index++, gp_Pnt(end->x(), end->y(), end->z()));  
+  Handle(Geom_BSplineCurve) NURBS = new Geom_BSplineCurve(ctrlPoints,_weights,_knots,_mult,degree,false);
+  TopoDS_Edge occEdge;
+  if (occv1 && occv2)
+    occEdge = BRepBuilderAPI_MakeEdge(NURBS,occv1->getShape(),occv2->getShape()).Edge();    
+  else
+    occEdge = BRepBuilderAPI_MakeEdge(NURBS).Edge();
+  return gm->_occ_internals->addEdgeToModel(gm,occEdge);
+  }
+  catch(Standard_Failure &err){
+    Msg::Error("%s", err.GetMessageString());
+  }
   return 0;
 }
 
-GEntity *OCCFactory::revolve(GModel *gm, GEntity* base, std::vector<double> p1, 
-                             std::vector<double> p2, double angle)
-{
+GEntity *OCCFactory::revolve (GModel *gm, GEntity* base,
+			      std::vector<double> p1, 
+			      std::vector<double> p2, double angle){
   if (!gm->_occ_internals)
     gm->_occ_internals = new OCC_Internals;
 
@@ -145,13 +225,11 @@ GEntity *OCCFactory::revolve(GModel *gm, GEntity* base, std::vector<double> p1,
     TopoDS_Solid result = TopoDS::Solid(MR.Shape());
     ret = gm->_occ_internals->addRegionToModel(gm, result);
   }
-  gm->glue(Precision::Confusion());
   return ret;
 }
 
-GEntity *OCCFactory::extrude(GModel *gm, GEntity* base, std::vector<double> p1, 
-                             std::vector<double> p2)
-{
+GEntity *OCCFactory::extrude (GModel *gm, GEntity* base,
+			      std::vector<double> p1, std::vector<double> p2){
   if (!gm->_occ_internals)
     gm->_occ_internals = new OCC_Internals;
 
@@ -182,7 +260,6 @@ GEntity *OCCFactory::extrude(GModel *gm, GEntity* base, std::vector<double> p1,
     TopoDS_Solid result = TopoDS::Solid(MP.Shape());
     ret = gm->_occ_internals->addRegionToModel(gm, result);    
   }
-  gm->glue(Precision::Confusion());
   return ret;
 }
 
@@ -192,11 +269,12 @@ GEntity *OCCFactory::addSphere(GModel *gm, double xc, double yc, double zc, doub
     gm->_occ_internals = new OCC_Internals;
 
   gp_Pnt aP(xc, yc, zc);  
-  gm->_occ_internals->buildShapeFromLists(BRepPrimAPI_MakeSphere(aP, radius).Shape());
+  TopoDS_Shape shape = BRepPrimAPI_MakeSphere(aP, radius).Shape();
+  gm->_occ_internals->buildShapeFromLists(shape);
   gm->destroy();
   gm->_occ_internals->buildLists();
   gm->_occ_internals->buildGModel(gm);  
-  return 0; // FIXME: should return GRegion?
+  return getOCCRegionByNativePtr(gm,TopoDS::Solid(shape));
 }
 
 GEntity *OCCFactory::addCylinder(GModel *gm, std::vector<double> p1,
@@ -223,11 +301,12 @@ GEntity *OCCFactory::addCylinder(GModel *gm, std::vector<double> p1,
     Msg::Error("Cylinder can't be computed from the given parameters");
     return 0;
   }
-  gm->_occ_internals->buildShapeFromLists(MC.Shape());
+  TopoDS_Shape shape = MC.Shape();
+  gm->_occ_internals->buildShapeFromLists(shape);
   gm->destroy();
   gm->_occ_internals->buildLists();
   gm->_occ_internals->buildGModel(gm);
-  return 0;  // FIXME: should return GRegion?
+  return getOCCRegionByNativePtr(gm,TopoDS::Solid(shape));
 }
 
 GEntity *OCCFactory::addTorus(GModel *gm, std::vector<double> p1, 
@@ -254,11 +333,12 @@ GEntity *OCCFactory::addTorus(GModel *gm, std::vector<double> p1,
     Msg::Error("Cylinder can't be computed from the given parameters");
     return 0;
   }
-  gm->_occ_internals->buildShapeFromLists(MC.Shape());
+  TopoDS_Shape shape = MC.Shape();
+  gm->_occ_internals->buildShapeFromLists(shape);
   gm->destroy();
   gm->_occ_internals->buildLists();
   gm->_occ_internals->buildGModel(gm);
-  return 0;  // FIXME: should return GRegion?
+  return getOCCRegionByNativePtr(gm,TopoDS::Solid(shape));
 }
 
 GEntity *OCCFactory::addCone(GModel *gm,  std::vector<double> p1, 
@@ -286,11 +366,12 @@ GEntity *OCCFactory::addCone(GModel *gm,  std::vector<double> p1,
     Msg::Error("Cylinder can't be computed from the given parameters");
     return 0;
   }
-  gm->_occ_internals->buildShapeFromLists(MC.Shape());
+  TopoDS_Shape shape = MC.Shape();
+  gm->_occ_internals->buildShapeFromLists(shape);
   gm->destroy();
   gm->_occ_internals->buildLists();
   gm->_occ_internals->buildGModel(gm);  
-  return 0; // FIXME: should return GRegion?
+  return getOCCRegionByNativePtr(gm,TopoDS::Solid(shape));
 }
 
 GEntity *OCCFactory::addBlock(GModel *gm, std::vector<double> p1, 
@@ -307,75 +388,145 @@ GEntity *OCCFactory::addBlock(GModel *gm, std::vector<double> p1,
     Msg::Error("Box can not be computed from the given point");
     return 0;
   }
-  gm->_occ_internals->buildShapeFromLists(MB.Shape());
+  TopoDS_Shape shape = MB.Shape();
+  gm->_occ_internals->buildShapeFromLists(shape);
   gm->destroy();
   gm->_occ_internals->buildLists();
   gm->_occ_internals->buildGModel(gm);
-  return 0;  // FIXME: should return GRegion?
+  return getOCCRegionByNativePtr(gm,TopoDS::Solid(shape));
 }
 
-GModel *OCCFactory::computeBooleanUnion(GModel *obj, GModel *tool, int createNewModel)
-{
-  OCC_Internals *occ_obj = obj->getOCCInternals();
-  OCC_Internals *occ_tool = tool->getOCCInternals();
-
-  if (!occ_obj || !occ_tool)return NULL;
-
-  if (createNewModel){
-    GModel *temp = new GModel;
-    temp->_occ_internals = new OCC_Internals;
-    temp->_occ_internals->addShapeToLists(occ_obj->getShape());
-    obj = temp;
+GModel *OCCFactory::computeBooleanUnion (GModel* obj, GModel* tool, int createNewModel){
+  try{ 
+    OCC_Internals *occ_obj  = obj->getOCCInternals();
+    OCC_Internals *occ_tool = tool->getOCCInternals();
+    
+    if (!occ_obj || !occ_tool)return NULL;
+    
+    if (createNewModel){
+      GModel *temp = new GModel;
+      temp->_occ_internals = new OCC_Internals;
+      temp->_occ_internals->addShapeToLists(occ_obj->getShape());
+      obj = temp;
+    }
+    obj->_occ_internals->applyBooleanOperator(occ_tool->getShape(),OCC_Internals::Fuse);
+    obj->destroy();
+    obj->_occ_internals->buildLists();
+    obj->_occ_internals->buildGModel(obj);
+  }
+  catch(Standard_Failure &err){
+    Msg::Error("%s", err.GetMessageString());
   }
-  obj->_occ_internals->applyBooleanOperator(occ_tool->getShape(), OCC_Internals::Fuse);
-  obj->destroy();
-  obj->_occ_internals->buildLists();
-  obj->_occ_internals->buildGModel(obj);
+    
   return obj;
 }
 
-GModel *OCCFactory::computeBooleanDifference(GModel *obj, GModel *tool,
-                                             int createNewModel)
-{
-  OCC_Internals *occ_obj = obj->getOCCInternals();
-  OCC_Internals *occ_tool = tool->getOCCInternals();
-
-  if (!occ_obj || !occ_tool)return NULL;
+GModel *OCCFactory::computeBooleanDifference (GModel* obj, GModel* tool, int createNewModel){
+  try{ 
+    OCC_Internals *occ_obj  = obj->getOCCInternals();
+    OCC_Internals *occ_tool = tool->getOCCInternals();
+    
+    if (!occ_obj || !occ_tool)return NULL;
+    
+    if (createNewModel){
+      GModel *temp = new GModel;
+      temp->_occ_internals = new OCC_Internals;
+      temp->_occ_internals->addShapeToLists(occ_obj->getShape());
+      obj = temp;
+    }
+    obj->getOCCInternals()->applyBooleanOperator(occ_tool->getShape(),OCC_Internals::Cut);
+    obj->destroy();
+    obj->_occ_internals->buildLists();
+    obj->_occ_internals->buildGModel(obj);
+  }
+  catch(Standard_Failure &err){
+    Msg::Error("%s", err.GetMessageString());
+  }
+  return obj;
+}
 
-  if (createNewModel){
-    GModel *temp = new GModel;
-    temp->_occ_internals = new OCC_Internals;
-    temp->_occ_internals->addShapeToLists(occ_obj->getShape());
-    obj = temp;
+GModel *OCCFactory::computeBooleanIntersection (GModel* obj, GModel* tool, int createNewModel){
+  try{
+    OCC_Internals *occ_obj  = obj->getOCCInternals();
+    OCC_Internals *occ_tool = tool->getOCCInternals();
+    
+    
+    if (!occ_obj || !occ_tool)return NULL;
+    
+    if (createNewModel){
+      GModel *temp = new GModel;
+      temp->_occ_internals = new OCC_Internals;
+      temp->_occ_internals->addShapeToLists(occ_obj->getShape());
+      obj = temp;
+    }
+    obj->getOCCInternals()->applyBooleanOperator(occ_tool->getShape(),OCC_Internals::Intersection);
+    obj->destroy();
+    obj->_occ_internals->buildLists();
+    obj->_occ_internals->buildGModel(obj);
+  }
+  catch(Standard_Failure &err){
+    Msg::Error("%s", err.GetMessageString());
   }
-  obj->getOCCInternals()->applyBooleanOperator(occ_tool->getShape(), 
-                                               OCC_Internals::Cut);
-  obj->destroy();
-  obj->_occ_internals->buildLists();
-  obj->_occ_internals->buildGModel(obj);
   return obj;
 }
 
-GModel *OCCFactory::computeBooleanIntersection(GModel *obj, GModel* tool, 
-                                               int createNewModel)
-{
-  OCC_Internals *occ_obj = obj->getOCCInternals();
-  OCC_Internals *occ_tool = tool->getOCCInternals();
+void OCCFactory::fillet (GModel *gm, std::vector<int> edges, double radius){
+  try{
+    std::vector<TopoDS_Edge> edgesToFillet;
+    for (int i=0;i<edges.size();i++){
+      GEdge *ed = gm->getEdgeByTag(edges[i]);
+      if (ed){
+	OCCEdge *occed = dynamic_cast<OCCEdge*>(ed);
+	if (occed)edgesToFillet.push_back(occed->getTopoDS_Edge());
+    }
+    }
+    gm->_occ_internals->Fillet(edgesToFillet,radius);
+    gm->destroy();
+    gm->_occ_internals->buildLists();
+    gm->_occ_internals->buildGModel(gm);  
+  }
+  catch(Standard_Failure &err){
+    Msg::Error("%s", err.GetMessageString());
+  }
+}
 
-  if (!occ_obj || !occ_tool)return NULL;
+void OCCFactory::translate (GModel *gm, std::vector<double> dx, int addToTheModel){
+  gp_Trsf transformation;
+  transformation.SetTranslation(gp_Pnt (0,0,0),gp_Pnt (dx[0],dx[1],dx[2]));
+  BRepBuilderAPI_Transform aTransformation(gm->_occ_internals->getShape(), transformation, Standard_False);
+  TopoDS_Shape temp = aTransformation.Shape();
+  if (!addToTheModel)gm->_occ_internals->loadShape(& temp);
+  else gm->_occ_internals->buildShapeFromLists(temp);
+  gm->destroy();
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);    
+}
 
-  if (createNewModel){
-    GModel *temp = new GModel;
-    temp->_occ_internals = new OCC_Internals;
-    temp->_occ_internals->addShapeToLists(occ_obj->getShape());
-    obj = temp;
-  }
-  obj->getOCCInternals()->applyBooleanOperator(occ_tool->getShape(),
-                                               OCC_Internals::Intersection);
-  obj->destroy();
-  obj->_occ_internals->buildLists();
-  obj->_occ_internals->buildGModel(obj);
-  return obj;
+void OCCFactory::rotate (GModel *gm, std::vector<double> p1,std::vector<double> p2, double angle, int addToTheModel){
+  const double x1 =p1[0]; 
+  const double y1 =p1[1]; 
+  const double z1 =p1[2]; 
+  const double x2 =p2[0]; 
+  const double y2 =p2[1]; 
+  const double z2 =p2[2]; 
+
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
+
+  gp_Trsf transformation;
+
+  gp_Vec direction (gp_Pnt (x1,y1,z1),gp_Pnt (x2,y2,z2));
+  gp_Ax1 axisOfRevolution (gp_Pnt (x1,y1,z1),direction);
+  transformation.SetRotation(axisOfRevolution, angle);
+  BRepBuilderAPI_Transform aTransformation(gm->_occ_internals->getShape(), transformation, Standard_False);
+  TopoDS_Shape temp = aTransformation.Shape();
+  if (!addToTheModel)gm->_occ_internals->loadShape(& temp);
+  else gm->_occ_internals->buildShapeFromLists(temp);
+  gm->destroy();
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);    
 }
 
+
 #endif
+
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index c4bbbc28ba..d58a6748e5 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -15,7 +15,6 @@
 class GVertex;
 class GEdge;
 class GModel;
-class binding;
 
 // Abstract CAD creation factory.
 class GModelFactory {
@@ -24,16 +23,22 @@ class GModelFactory {
   virtual ~GModelFactory(){}
   // brep primitives
   enum arcCreationMethod {THREE_POINTS=1, CENTER_START_END=2};
-  enum splineType {BEZIER=1, CATMULL_ROM=2};
+  enum splineType {BEZIER=1, BSPLINE=2};
   virtual GVertex *addVertex(GModel *gm, double x, double y, double z, 
                              double lc) = 0;
   virtual GEdge *addLine(GModel *, GVertex *v1, GVertex *v2) = 0;
   virtual GEdge *addCircleArc(GModel *gm, const arcCreationMethod &method,
                               GVertex *start, GVertex *end, 
                               const SPoint3 &aPoint) = 0;
-  virtual GEdge *addSpline(GModel *gm, const splineType &type,
-                           GVertex *start, GVertex *end, 
-                           fullMatrix<double> *controlPoints) = 0;
+  virtual GEdge *addSpline (GModel *gm,const splineType &type,
+			    GVertex *start, 
+			    GVertex *end, 
+			    std::vector<std::vector<double> > controlPoints) = 0;
+  virtual GEdge *addNURBS(GModel *gm, GVertex *start, GVertex *end,
+			  std::vector<std::vector<double> > controlPoints, 
+			  std::vector<double> knots,
+			  std::vector<double> weights, 
+			  std::vector<int> multiplicity) = 0;
   virtual GEntity *revolve(GModel *gm, GEntity*, std::vector<double> p1, 
                            std::vector<double> p2, double angle) = 0;
   virtual GEntity *extrude(GModel *gm, GEntity*, std::vector<double> p1, 
@@ -53,52 +58,61 @@ class GModelFactory {
                            std::vector<double> p2, double radius1,
                            double radius2) = 0; 
 
-  // boolean operators acting on two GModels
-  virtual GModel *computeBooleanUnion(GModel *obj, GModel *tool, 
-                                      int createNewModel) = 0;
-  virtual GModel *computeBooleanIntersection(GModel *obj, GModel *tool,
-                                             int createNewModel) = 0;
-  virtual GModel *computeBooleanDifference(GModel *obj, GModel *tool,
-                                           int createNewModel) = 0;
+  // here, we should give a list of GEdges. Yet, I still can't figure out how
+  // to get those automatically ... wait and see
+  virtual void fillet (GModel *gm, std::vector<int> edges, double radius) = 0; 
+
+  // rigid body motions
+  virtual void translate (GModel *gm, std::vector<double> dx, int addToTheModel) = 0;
+  virtual void rotate (GModel *gm, std::vector<double> p1,std::vector<double> p2, double angle, int addToTheModel) = 0;
+
+  // boolean operators acting on 2 GEntities
+  virtual GModel * computeBooleanUnion (GModel *obj, GModel*tool, int createNewModel) = 0;
+  virtual GModel * computeBooleanIntersection (GModel *obj, GModel*tool, int createNewModel) = 0;
+  virtual GModel * computeBooleanDifference (GModel *obj, GModel*tool, int createNewModel) = 0;
 };
 
 #if defined(HAVE_OCC)
 
 class OCCFactory : public GModelFactory {
  public:
-  OCCFactory(){}
-  virtual ~OCCFactory(){}
-  virtual GVertex *addVertex(GModel *gm, double x, double y, double z, 
-                             double lc);
-  virtual GEdge *addLine(GModel *, GVertex *v1, GVertex *v2);
-  virtual GEdge *addCircleArc(GModel *gm, const arcCreationMethod &method,
-                              GVertex *start, GVertex *end, 
-                              const SPoint3 &aPoint);
-  virtual GEdge *addSpline(GModel *gm, const splineType &type,
-                           GVertex *start, GVertex *end, 
-                           fullMatrix<double> *controlPoints);
-  virtual GEntity *revolve(GModel *gm, GEntity*, std::vector<double> p1, 
-                           std::vector<double> p2, double angle);
-  virtual GEntity *extrude(GModel *gm, GEntity*, std::vector<double> p1, 
-                           std::vector<double> p2);
-  virtual GEntity *addSphere(GModel *gm, double cx, double cy, double cz, 
-                             double radius) ; 
-  virtual GEntity *addCylinder(GModel *gm, std::vector<double> p1, 
-                               std::vector<double> p2, double radius);
-  virtual GEntity *addTorus(GModel *gm, std::vector<double> p1, 
-                            std::vector<double> p2, double radius1, 
-                            double radius2);
-  virtual GEntity *addBlock(GModel *gm, std::vector<double> p1,
-                            std::vector<double> p2);
-  virtual GEntity *addCone(GModel *gm, std::vector<double> p1,
-                           std::vector<double> p2, double radius1,
-                           double radius2); 
-  virtual GModel *computeBooleanUnion(GModel *obj, GModel *tool, 
-                                      int createNewModel);
-  virtual GModel *computeBooleanIntersection(GModel *obj, GModel *tool,
-                                             int createNewModel);
-  virtual GModel *computeBooleanDifference(GModel *obj, GModel *tool,
-                                           int createNewModel);
+  OCCFactory (){}
+  GVertex *addVertex (GModel *gm,double x, double y, double z, double lc);
+  virtual GEdge *addLine (GModel *gm,GVertex *v1, GVertex *v2);
+  GEdge *addCircleArc (GModel *gm,const arcCreationMethod &method,				  
+		       GVertex *start, 
+		       GVertex *end, 
+		       const SPoint3 &aPoint);
+  GEdge *addSpline (GModel *gm,const splineType &type,
+		    GVertex *start, 
+		    GVertex *end, 
+		    std::vector<std::vector<double> > controlPoints);
+  GEdge *addNURBS(GModel *gm,
+		  GVertex *start, GVertex *end,
+		  std::vector<std::vector<double> > controlPoints, 
+		  std::vector<double> knots,
+		  std::vector<double> weights, 
+		  std::vector<int> multiplicity);
+  
+  GEntity* revolve (GModel *gm, GEntity*,std::vector<double> p1,std::vector<double> p2, double angle);
+  GEntity* extrude (GModel *gm, GEntity*,std::vector<double> p1,std::vector<double> p2);
+
+  GEntity * addSphere   (GModel *gm,double cx, double cy, double cz, double radius); 
+  GEntity * addCylinder (GModel *gm,std::vector<double> p1, std::vector<double> p2, 
+		      double radius); 
+  GEntity * addTorus   (GModel *gm,std::vector<double> p1, std::vector<double> p2, 
+		     double radius1, double radius2); 
+  GEntity * addBlock  (GModel *gm,std::vector<double> p1, std::vector<double> p2); 
+  GEntity * addCone   (GModel *gm,std::vector<double> p1, std::vector<double> p2, double radius1, double radius2); 
+  // rigid body motions
+  void translate (GModel *gm, std::vector<double> dx, int addToTheModel);
+  void rotate (GModel *gm, std::vector<double> p1,std::vector<double> p2, double angle, int addToTheModel);
+  // booleans
+  GModel * computeBooleanUnion (GModel *obj, GModel*tool, int createNewModel);
+  GModel * computeBooleanIntersection (GModel *obj, GModel*tool, int createNewModel);
+  GModel * computeBooleanDifference (GModel *obj, GModel*tool, int createNewModel);
+  // advanced
+  void fillet (GModel *gm, std::vector<int> edges, double radius); 
 };
 
 #endif
diff --git a/Geo/GModelIO_ACIS.cpp b/Geo/GModelIO_ACIS.cpp
new file mode 100644
index 0000000000..7052d0e7d5
--- /dev/null
+++ b/Geo/GModelIO_ACIS.cpp
@@ -0,0 +1,149 @@
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "GmshMessage.h"
+#include "GModelIO_ACIS.h"
+#include "ACISVertex.h"
+#include "ACISEdge.h"
+
+#if defined(HAVE_ACIS)
+#define MacX 1
+#define mac 1
+#define ANSI 1
+
+#include <acis.hxx>
+#include <base.hxx>
+#include <license.hxx>
+#include <spa_unlock_result.hxx>
+#include <spa_lic_err_gui.hxx>
+#include <api.hxx>
+#include <kernapi.hxx>
+#include <mmgr.err>
+#include <lists.hxx>
+#include <acistype.hxx>
+
+
+class ACIS_Internals {
+public:
+  ENTITY_LIST entities;
+  ACIS_Internals();
+  ~ACIS_Internals();
+  void loadSAT ( std::string fileName, GModel*);  
+  void addVertices (GModel *gm, ENTITY_LIST &l);
+  void addEdges (GModel *gm, ENTITY_LIST &l);
+  //  void addFaces (GModel *gm, ENTITY_LIST &l);
+};
+
+ACIS_Internals::ACIS_Internals() {
+  // put your acis unlock string here ...
+#include "ACISLICENSE.h"
+
+  spa_unlock_result out = spa_unlock_products( unlock_str );
+
+  outcome prout = api_start_modeller(0);
+  if (!prout.ok()){
+    Msg::Error("Unable to start ACIS");
+  }
+}
+ACIS_Internals::~ACIS_Internals(){
+  outcome prout = api_stop_modeller();
+  if (!prout.ok()){
+    Msg::Error("Unable to stop ACIS");
+  }  
+}
+
+void ACIS_Internals :: addVertices (GModel *gm, ENTITY_LIST &l){
+  l.init();
+  ENTITY *e;
+  while(e = l.next()){
+    VERTEX *av = dynamic_cast<VERTEX*>(e);
+    if (av){
+      GVertex *v = getACISVertexByNativePtr(gm,av);
+      if (!v)
+	gm->add(new ACISVertex (gm,gm->maxVertexNum()+1,av));
+    }
+  }
+}
+
+void ACIS_Internals :: addEdges (GModel *gm, ENTITY_LIST &l){
+  l.init();
+  ENTITY *e;
+  while(e = l.next()){
+    EDGE *av = dynamic_cast<EDGE*>(e);
+    if (av){
+      GEdge *v = getACISEdgeByNativePtr(gm,av);
+      if (!v){
+	GVertex *v1 = getACISVertexByNativePtr(gm,av->start());
+	GVertex *v2 = getACISVertexByNativePtr(gm,av->end());      
+	gm->add(new ACISEdge(gm,av,gm->maxEdgeNum()+1,v1,v2));
+      }
+    }
+  }
+}
+/*
+void ACIS_Internals :: addFaces (GModel *gm, ENTITY_LIST &l){
+  l.init();
+  ENTITY *e;
+  while(e = l.next()){
+    FACE *av = dynamic_cast<FACE*>(e);
+    if (av){
+      GFace *v = getACISFaceByNativePtr(gm,av);
+      if (!v){
+	GVertex *v1 = getACISVertexByNativePtr(gm,av->start());
+	GVertex *v2 = getACISVertexByNativePtr(gm,av->end());      
+	gm->add(new ACISEdge(gm,av,gm->maxEdgeNum()+1,v1,v2));
+      }
+    }
+  }
+}
+*/
+
+void ACIS_Internals::loadSAT (std::string fileName, GModel *gm) {
+  FILE *f = fopen (fileName.c_str(), "r");
+  if (!f){
+    return;
+  }  
+  outcome prout = api_restore_entity_list(f,1,entities); 	
+  if (!prout.ok()){
+    Msg::Error("Unable to load ACIS FILE %d",fileName.c_str());
+    return;
+  }  
+  Msg::Info("ACIS FILE %d Loaded",fileName.c_str());
+  
+  ENTITY *e;
+  entities.init();
+  while(e = entities.next()){
+    printf("an entity\n");
+    if (is_VERTEX(e)){
+      printf("VERTEX FOUND\n");
+    }
+    if (is_BODY(e)){
+      {
+	ENTITY_LIST vertex_list;
+	outcome prout = api_get_vertices (e,vertex_list);
+	addVertices (gm,vertex_list);
+	printf("BODY COUNT %d !\n",vertex_list.count());
+      }
+      {
+	ENTITY_LIST edge_list;
+	outcome prout = api_get_edges (e,edge_list);
+	addEdges (gm,edge_list);
+	printf("BODY COUNT %d !\n",edge_list.count());
+      }
+    }
+  }
+}
+
+int GModel::readACISSAT(const std::string &fn)
+{
+  _acis_internals = new ACIS_Internals;
+  _acis_internals->loadSAT(fn,this);
+  return 1;
+}
+#else
+int GModel::readACISSAT(const std::string &fn)
+{
+  Msg::Error("Gmsh must be compiled with ACIS support to load '%s'",
+             fn.c_str());
+  return 0;
+}
+#endif
diff --git a/Geo/GModelIO_ACIS.h b/Geo/GModelIO_ACIS.h
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 4ff8a85ad5..fef6327e6c 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -18,6 +18,8 @@
 
 #if defined(HAVE_OCC)
 
+void addSimpleShapes(TopoDS_Shape theShape, TopTools_ListOfShape &theList);
+
 void OCC_Internals::buildLists()
 {
   somap.Clear();
@@ -34,7 +36,13 @@ void OCC_Internals::buildShapeFromLists(TopoDS_Shape _shape)
   BRep_Builder B;
   TopoDS_Compound C;
   B.MakeCompound(C);
-  B.Add(C,_shape);
+  //  B.Add(C,_shape);
+
+  TopTools_ListOfShape theList;
+  addSimpleShapes(_shape, theList);
+  TopTools_ListIteratorOfListOfShape itSub1 (theList);
+  for (; itSub1.More(); itSub1.Next()) B.Add(C,itSub1.Value());
+
   for(int i = 1; i <= vmap.Extent(); i++) B.Add(C, vmap(i));
   for(int i = 1; i <= emap.Extent(); i++) B.Add(C, emap(i));
   for(int i = 1; i <= fmap.Extent(); i++) B.Add(C, fmap(i));
@@ -463,107 +471,66 @@ GVertex *OCC_Internals::addVertexToModel(GModel *model, TopoDS_Vertex vertex)
 {
   GVertex *gv = getOCCVertexByNativePtr(model, vertex);
   if (gv) return gv;
-  buildShapeFromLists(vertex);
-  gv = new OCCVertex(model, model->maxVertexNum() + 1, vertex);
-  model->add(gv);
-  return gv;
+  addShapeToLists(vertex);
+  buildGModel(model);
+  return getOCCVertexByNativePtr (model,vertex);
 }
 
 GEdge *OCC_Internals::addEdgeToModel(GModel *model, TopoDS_Edge edge)
 {
   GEdge *ge = getOCCEdgeByNativePtr(model, edge);
   if (ge) return ge;
-  buildShapeFromLists(edge);
-  TopoDS_Vertex occv1 = TopExp::FirstVertex(edge);
-  TopoDS_Vertex occv2 = TopExp::LastVertex(edge);
-  GVertex *v1 = addVertexToModel(model, occv1);
-  GVertex *v2 = addVertexToModel(model, occv2);
-  OCCEdge *e = new OCCEdge(model, edge, model->maxEdgeNum() + 1, v1, v2);
-  model->add(e);
-  return e;
-}
-
-GEdge *OCC_Internals::addEdgeToModel(GModel *model, TopoDS_Edge edge, GVertex *g1, 
-                                     GVertex *g2)
-{
-  OCCEdge *e = new OCCEdge(model, edge, model->maxEdgeNum() + 1, g1, g2);
-  e->replaceEndingPoints (g1,g2);
-  model->add(e);
-  return e;
+  addShapeToLists(edge);
+  buildGModel(model);
+  return getOCCEdgeByNativePtr(model,edge);
 }
 
-GFace *OCC_Internals::addFaceToModel(GModel *model, TopoDS_Face face, int i)
-{
+GFace* OCC_Internals::addFaceToModel(GModel *model, TopoDS_Face face){
   GFace *gf = getOCCFaceByNativePtr(model,face);
   if (gf) return gf;
-
-  if (i < 0){
-    std::list<GEdge*> _edges;
-    TopExp_Explorer exp2, exp3;
-    for(exp2.Init(face, TopAbs_WIRE); exp2.More(); exp2.Next()){
-      TopoDS_Wire wire = TopoDS::Wire(exp2.Current());
-      for(exp3.Init(wire, TopAbs_EDGE); exp3.More(); exp3.Next()){          
-	TopoDS_Edge edge = TopoDS::Edge(exp3.Current());
-	_edges.push_back(addEdgeToModel(model, edge));
-      }
-    }
-    i = model->maxFaceNum() + 1;
-    OCCFace *f = new OCCFace(model, face, i);
-    model->add(f);
-    model->glue(Precision::Confusion());
-    return f;
-  }
-  OCCFace *f = new OCCFace(model, face, i);
-  model->add(f);
-  return f;
+  addShapeToLists(face);
+  buildGModel(model);
+  return getOCCFaceByNativePtr(model,face);
 }
 
-GRegion *OCC_Internals::addRegionToModel(GModel *model, TopoDS_Solid region, int i)
-{
-  if (i < 0){
-    TopExp_Explorer exp0, exp1, exp2;
-    for(exp0.Init(region, TopAbs_SOLID); exp0.More(); exp0.Next()){
-      TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
-      for(exp1.Init(exp0.Current(), TopAbs_SHELL); exp1.More(); exp1.Next()){
-	TopoDS_Shell shell = TopoDS::Shell(exp1.Current());
-	for(exp2.Init(shell, TopAbs_FACE); exp2.More(); exp2.Next()){
-	  TopoDS_Face face = TopoDS::Face(exp2.Current());
-	  addFaceToModel(model, face, -1);
-	}
-      }
-    }
-    i = model->maxRegionNum() + 1;
-  }
-  OCCRegion *r = new OCCRegion(model, region, i);
-  model->add(r);
-  return r;
+GRegion* OCC_Internals::addRegionToModel(GModel *model, TopoDS_Solid region){
+  
+  GRegion *gr  = getOCCRegionByNativePtr(model,region);
+  if (gr) return gr;
+  addShapeToLists(region);
+  //  buildLists();
+  // TEST
+  buildGModel(model);
+  return getOCCRegionByNativePtr(model,region);
 }
 
 void OCC_Internals::buildGModel(GModel *model)
 {
   // building geom vertices
   int nvertices = vmap.Extent();
-  printf("%d vertices\n",nvertices);
   for(int i = 1; i <= nvertices; i++){
-    model->add(new OCCVertex(model, i, TopoDS::Vertex(vmap(i))));
+    if (!getOCCVertexByNativePtr(model,TopoDS::Vertex(vmap(i))))
+      model->add(new OCCVertex(model, i, TopoDS::Vertex(vmap(i))));
   }
   // building geom edges
   int nedges = emap.Extent();
   for(int i = 1; i <= nedges; i++){
     int i1 = vmap.FindIndex(TopExp::FirstVertex(TopoDS::Edge(emap(i)))); 
     int i2 = vmap.FindIndex(TopExp::LastVertex(TopoDS::Edge(emap(i))));
-    model->add(new OCCEdge(model, TopoDS::Edge(emap(i)), i,
-                           model->getVertexByTag(i1), model->getVertexByTag(i2)));
+    if (!getOCCEdgeByNativePtr(model,TopoDS::Edge(emap(i))))
+      model->add(new OCCEdge(model, TopoDS::Edge(emap(i)), i, model->getVertexByTag(i1), model->getVertexByTag(i2)));
   }
   // building geom faces
   int nfaces = fmap.Extent();
   for(int i = 1; i <= nfaces; i++){
-    model->add(new OCCFace(model, TopoDS::Face(fmap(i)), i));
+    if (!getOCCFaceByNativePtr(model,TopoDS::Face(fmap(i))))
+      model->add(new OCCFace(model, TopoDS::Face(fmap(i)), i));
   }
   // building geom regions
   int nvolumes = somap.Extent();
   for(int i = 1; i <= nvolumes; i++){
-    model->add(new OCCRegion(model, TopoDS::Solid(somap(i)), i));
+    if (!getOCCRegionByNativePtr(model,TopoDS::Solid(somap(i))))
+      model->add(new OCCRegion(model, TopoDS::Solid(somap(i)), i));
   }
 }
 
@@ -592,53 +559,6 @@ void addSimpleShapes(TopoDS_Shape theShape, TopTools_ListOfShape &theList)
   }
 }
 
-static TopoDS_Shape GlueFaces(const TopoDS_Shape &theShape,
-                              const Standard_Real theTolerance)
-{
-  Msg::Error("glue !");
-  return theShape;
-
-//   Standard_Integer iErr, iWrn;
-//   TopoDS_Shape aRes;
-//   GEOMAlgo_Gluer aGluer;
-
-//   aGluer.SetShape(theShape);
-//   aGluer.SetTolerance(theTolerance);
-//   aGluer.SetCheckGeometry(Standard_True);
-
-//   aGluer.Perform();
-
-//   iErr = aGluer.ErrorStatus();
-//   if (iErr) {
-//     switch (iErr) {
-//     case 2:
-//       Msg::Error("No vertices found in source shape");
-//       break;
-//     case 5:
-//       Msg::Error("Source shape is Null");
-//       break;
-//     case 6:
-//       Msg::Error("Result shape is Null");
-//       break;
-//     case 200:
-//       Msg::Error("Error occured during check of geometric coincidence");
-//       break;
-//     default:
-//       {
-//         // description of all errors see in GEOMAlgo_Gluer.cxx
-//         TCollection_AsciiString aMsg ("Error in GEOMAlgo_Gluer with code ");
-//         aMsg += TCollection_AsciiString(iErr);
-//         Msg::Error(aMsg.ToCString());
-//         break;
-//       }
-//     }
-//     return aRes;
-//   }
-
-//   aRes = aGluer.Result();
-
-//   return aRes;
-}
 
 void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool, const BooleanOperator &op)
 {
@@ -693,9 +613,9 @@ void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool, const BooleanOperato
             TopoDS_Shape aValueC = itSubC.Value();
             if (aValueC.ShapeType() != TopAbs_SOLID) isOnlySolids = false;
           }
-          if (isOnlySolids)
-            theNewShape = GlueFaces(C, Precision::Confusion());
-          else
+	  //          if (isOnlySolids)
+	  //            theNewShape = GlueFaces(C, Precision::Confusion());
+	  //          else
             theNewShape = C;
         }
         shape = theNewShape;
@@ -752,10 +672,10 @@ void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool, const BooleanOperato
             TopoDS_Shape aValueC = itSubC.Value();
             if (aValueC.ShapeType() != TopAbs_SOLID) isOnlySolids = false;
           }
-          if (isOnlySolids)
-            theNewShape = GlueFaces(C, Precision::Confusion());
-          else
-            theNewShape = C;
+	  //          if (isOnlySolids)
+	  //            theNewShape = GlueFaces(C, Precision::Confusion());
+	  //          else
+	  theNewShape = C;
         }
         shape = theNewShape;
       }      
@@ -823,7 +743,7 @@ void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool, const BooleanOperato
   }
 }
   
-void OCC_Internals::fillet(std::vector<TopoDS_Edge> &edgesToFillet,
+void OCC_Internals::Fillet(std::vector<TopoDS_Edge> &edgesToFillet,
                            double Radius)
 {
   // create a tool for fillet
diff --git a/Geo/GModelIO_OCC.h b/Geo/GModelIO_OCC.h
index 168360d4d0..d986ec8a94 100644
--- a/Geo/GModelIO_OCC.h
+++ b/Geo/GModelIO_OCC.h
@@ -42,13 +42,23 @@ class OCC_Internals {
   void writeSTEP(const char *);  
   void loadIGES(const char *);
   void loadShape(const TopoDS_Shape *);
-  GVertex *addVertexToModel(GModel *model, TopoDS_Vertex v);
-  GEdge *addEdgeToModel(GModel *model, TopoDS_Edge e);
-  GEdge *addEdgeToModel(GModel *model, TopoDS_Edge e, GVertex *v1, GVertex *v2);
-  GFace *addFaceToModel(GModel *model, TopoDS_Face f, int i=-1);
-  GRegion *addRegionToModel(GModel *model, TopoDS_Solid r, int i=-1);
   void buildGModel(GModel *gm);
-  void fillet(std::vector<TopoDS_Edge> &shapes, double radius);
+  GVertex * addVertexToModel(GModel *model, TopoDS_Vertex v);
+  GEdge   * addEdgeToModel  (GModel *model, TopoDS_Edge e);
+  GFace   * addFaceToModel  (GModel *model, TopoDS_Face f);
+  GRegion * addRegionToModel (GModel *model, TopoDS_Solid r);
+
+  void Box(const SPoint3 &p1, const SPoint3 &p2, const BooleanOperator &op);
+  void Sphere(const SPoint3 &center, const double &radius, const BooleanOperator &op);
+  void Cylinder(const SPoint3 &bottom_center, const SVector3 &dir, double R, double H,
+                const BooleanOperator &op);
+  void Cone(const SPoint3 &bottom_center, const SVector3 &dir, double R1, double R2, 
+             double H,  const BooleanOperator &op);
+  void Torus(const SPoint3 &bottom_center, const SVector3 &dir, double R1, double R2, 
+             const BooleanOperator &op);
+  void Torus(const SPoint3 &bottom_center, const SVector3 &dir, double R1, double R2, 
+             double angle,  const BooleanOperator &op);
+  void Fillet(std::vector<TopoDS_Edge> &shapes, double radius);
   void applyBooleanOperator(TopoDS_Shape tool, const BooleanOperator &op);
 };
 
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index b5320c0c32..3e6c9b6a2c 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -99,4 +99,11 @@ void GVertex::registerBindings(binding *b)
   classBinding *cb = b->addClass<GVertex>("GVertex");
   cb->setDescription("A GVertex is a geometrical 0D entity");
   cb->setParentClass<GEntity>();
+  methodBinding *cm;
+  cm = cb->addMethod("x", (double (GVertex::*)() const) &GVertex::x);
+  cm->setDescription("Return the x-coordinate.");
+  cm = cb->addMethod("y", (double (GVertex::*)() const) &GVertex::y);
+  cm->setDescription("Return the y-coordinate.");
+  cm = cb->addMethod("z", (double (GVertex::*)() const) &GVertex::z);
+  cm->setDescription("Return the z-coordinate.");
 }
diff --git a/Geo/OCCIncludes.h b/Geo/OCCIncludes.h
index 80680b4b30..e6f54ff398 100644
--- a/Geo/OCCIncludes.h
+++ b/Geo/OCCIncludes.h
@@ -65,7 +65,6 @@ using std::iostream;
 #include "BRepAdaptor_Curve.hxx"
 #include "Poly_Triangulation.hxx"
 #include "Poly_Array1OfTriangle.hxx"
-#include "TColgp_Array1OfPnt2d.hxx"
 #include "Poly_Triangle.hxx"
 #include "GProp_GProps.hxx"
 #include "BRepGProp.hxx"
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index 305452081e..946da65ffe 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -114,4 +114,17 @@ void OCCRegion::replaceFacesInternal(std::list<GFace*> &new_faces)
   setup();
 }
 
+GRegion *getOCCRegionByNativePtr(GModel *model, TopoDS_Solid toFind)
+{
+  GModel::riter it =model->firstRegion();
+  for (; it !=model->lastRegion(); it++){
+    OCCRegion *occr = dynamic_cast<OCCRegion*>(*it);
+    if (occr){
+      if (toFind.IsSame(occr->getTopoDS_Shape())){
+	return *it;
+      }
+    }
+  }
+  return 0;
+}
 #endif
diff --git a/Geo/OCCRegion.h b/Geo/OCCRegion.h
index 33b5a1de18..60cdac47e9 100644
--- a/Geo/OCCRegion.h
+++ b/Geo/OCCRegion.h
@@ -25,6 +25,7 @@ class OCCRegion : public GRegion {
   TopoDS_Solid getTopoDS_Shape() {return s;}
 };
 
+GRegion *getOCCRegionByNativePtr(GModel *model, TopoDS_Solid toFind);
 #endif
 
 #endif
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index fa470e20ef..8ccca3e8f4 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -62,7 +62,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
     // database (if any)
     _try(MatSetFromOptions(_a));
     // preallocation option must be set after other options
-    PetscInt prealloc = 200;
+    PetscInt prealloc = 300;
     PetscTruth set;
     PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
     _try(MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL));
@@ -146,7 +146,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
     // set some default options
     _try(PCSetType(pc, PCILU));
     _try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
-    _try(PCFactorSetLevels(pc, 6));
+    _try(PCFactorSetLevels(pc, 1));
     _try(KSPSetTolerances(ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
     // override the default options with the ones from the option
     // database (if any)
-- 
GitLab