From 945ff17bb5b27543ab2a41b7efb63d18105e7878 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 29 Apr 2010 06:47:20 +0000
Subject: [PATCH] merged JF and my new solid modeler code

---
 Common/LuaBindings.cpp           |  26 +-
 Geo/GModel.cpp                   | 337 +++++++++++++++++++-------
 Geo/GModel.h                     |  65 +++--
 Geo/GModelFactory.cpp            | 400 +++++++++++++------------------
 Geo/GModelFactory.h              | 183 ++++++--------
 Geo/GModelIO_OCC.cpp             |  34 ++-
 Geo/GModelIO_OCC.h               |   2 +-
 Solver/linearSystemPETSc.cpp     |   5 +-
 benchmarks/boolean/square1.lua   |  11 +
 benchmarks/boolean/square2.lua   |  15 ++
 benchmarks/boolean/wikipedia.lua |  26 ++
 11 files changed, 634 insertions(+), 470 deletions(-)
 create mode 100644 benchmarks/boolean/square1.lua
 create mode 100644 benchmarks/boolean/square2.lua
 create mode 100644 benchmarks/boolean/wikipedia.lua

diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 5ceb187b06..4f1ae0e16a 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -26,16 +26,15 @@
 #include "drawContext.h"
 #include "GmshMessage.h"
 #include "linearSystem.h"
-#include "GModelFactory.h"
 
 #if defined(HAVE_SOLVER)
 #include "linearSystemCSR.h"
 #endif
 
 extern "C" {
-  #include "lua.h"
-  #include "lualib.h"
-  #include "lauxlib.h"
+#include "lua.h"
+#include "lualib.h"
+#include "lauxlib.h"
 }
 
 #if defined(HAVE_READLINE)
@@ -49,7 +48,7 @@ class gmshOptions {
   gmshOptions(){}
   void colorSet(std::string category, int index, std::string name, int value)
   {
-    GmshSetOption(category,name,(unsigned int)(value), index);
+    GmshSetOption(category, name, (unsigned int)(value), index);
   }
   int colorGet(std::string category, int index, std::string name)
   {
@@ -77,7 +76,8 @@ class gmshOptions {
   {
     GmshSetOption(category, name, value, index);
   }
-  static void registerBindings(binding *b) {
+  static void registerBindings(binding *b)
+  {
     classBinding *cb = b->addClass<gmshOptions>("gmshOptions");
     cb->setDescription("access the gmsh option database");
     methodBinding *mb;
@@ -227,7 +227,7 @@ static int luaClear(lua_State *L)
 }
 #endif
 
-static int luaRefresh (lua_State *L)
+static int luaRefresh(lua_State *L)
 {
   drawContext::global()->draw();
   return 0;
@@ -371,11 +371,11 @@ binding::binding()
   // Register Lua bindings
   DocRecord::registerBindings(this);
   GEntity::registerBindings(this);
+  GVertex::registerBindings(this);
   GEdge::registerBindings(this);
   GFace::registerBindings(this);
-  GModel::registerBindings(this);
   GRegion::registerBindings(this);
-  GVertex::registerBindings(this);
+  GModel::registerBindings(this);
   MElement::registerBindings(this);
   MVertex::registerBindings(this);
   fullMatrix<double>::registerBindings(this);
@@ -386,12 +386,14 @@ binding::binding()
   function::registerBindings(this);
   linearSystemCSRGmm<double>::registerBindings(this);
 #endif
-  GModelFactory::registerBindings(this);
 }
 
 binding *binding::_instance=NULL;
-binding::~binding() {
-  for (std::map<std::string,classBinding *>::iterator it =  classes.begin(); it != classes.end(); it++) {
+
+binding::~binding()
+{
+  for (std::map<std::string,classBinding *>::iterator it = classes.begin(); 
+       it != classes.end(); it++) {
     delete it->second;
   }
 }
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 7332f42954..4f953e58b3 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -8,6 +8,7 @@
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "GModel.h"
+#include "GModelFactory.h"
 #include "MPoint.h"
 #include "MLine.h"
 #include "MTriangle.h"
@@ -41,13 +42,16 @@ int GModel::_current = -1;
 
 GModel::GModel(std::string name)
   : _name(name), _visible(1), _octree(0),
-    _geo_internals(0), _occ_internals(0), _fm_internals(0),
+    _geo_internals(0), _occ_internals(0), _fm_internals(0), _factory(0),
     _fields(0), _currentMeshEntity(0), normals(0)
 {
   partitionSize[0] = 0; partitionSize[1] = 0;
   list.push_back(this);
   // at the moment we always create (at least an empty) GEO model
   _createGEOInternals();
+#if defined(HAVE_OCC)
+  _factory = new OCCFactory();
+#endif
 #if defined(HAVE_MESH)
   _fields = new FieldManager();
 #endif
@@ -84,7 +88,23 @@ int GModel::setCurrent(GModel *m)
       break;
     }
   }
-        return _current;
+  return _current;
+}
+
+void GModel::setFactory(std::string name)
+{
+  if(_factory) delete _factory;
+  _factory = 0;
+  if(name == "Gmsh") {
+    Msg::Error("Gmsh factory not created yet");
+  }
+  else if(name == "OpenCASCADE"){
+#if defined(HAVE_OCC)
+    _factory = new OCCFactory();
+#else
+    Msg::Error("Missing OpenCASCADE support");
+#endif
+  }
 }
 
 GModel *GModel::findByName(std::string name)
@@ -1402,18 +1422,124 @@ void GModel::save(std::string fileName)
   GModel::setCurrent(temp);
 }
 
-static void ComputeDuplicates(GModel * model, 
-                              std::multimap<GVertex*,GVertex*>  & Unique2Duplicates,
-                              std::map<GVertex*,GVertex*> & Duplicates2Unique,
-                              const double &eps)
+GVertex *GModel::addVertex(double x, double y, double z, double lc)
+{
+  if(_factory) return _factory->addVertex(this, x, y, z, lc);
+  return 0;
+}
+
+GEdge *GModel::addLine(GVertex *v1, GVertex *v2)
+{
+  if(_factory) return _factory->addLine(this, v1, v2);
+  return 0;
+}
+
+GEdge *GModel::addCircleArcCenter(double x, double y, double z, GVertex *start, 
+                                  GVertex *end)
+{
+  if(_factory)
+    return _factory->addCircleArc(this, GModelFactory::CENTER_START_END,
+                                  start, end, SPoint3(x, y, z));
+  return 0;
+}
+
+GEdge *GModel::addCircleArc3Points(double x, double y, double z, GVertex *start,
+                                   GVertex *end)
+{
+  if(_factory) 
+    return _factory->addCircleArc(this, GModelFactory::THREE_POINTS, 
+                                  start, end, SPoint3(x, y, z));
+  return 0;
+}
+
+GEdge *GModel::addBezier(GVertex *start, GVertex *end, 
+                         fullMatrix<double> *controlPoints)
+{
+  if(_factory) 
+    return _factory->addSpline(this, GModelFactory::BEZIER, start, end, 
+                               controlPoints);
+  return 0;
+}
+
+GEntity *GModel::revolve(GEntity *e, double angle, fullMatrix<double> *axis)
+{
+  if(_factory)
+    return _factory->revolve(this, e, (*axis)(0, 0), (*axis)(0, 1), (*axis)(0, 2),
+                             (*axis)(1, 0), (*axis)(1, 1), (*axis)(1, 2), angle);
+  return 0;
+}
+
+GEntity *GModel::extrude(GEntity *e, fullMatrix<double> *axis)
+{
+  if(_factory) 
+    return _factory->extrude(this, e, (*axis)(0, 0), (*axis)(0, 1), (*axis)(0, 2),
+                             (*axis)(1, 0), (*axis)(1, 1), (*axis)(1, 2));
+  return 0;
+}
+
+GEntity *GModel::addSphere(double cx, double cy, double cz, double radius)
+{
+  if(_factory) return _factory->addSphere(this, cx, cy, cz, radius);
+  return 0;
+}
+
+GEntity *GModel::addCylinder(std::vector<double> p1, std::vector<double> p2, 
+                             double radius)
 {
-  // in a first time, we use a greedy algorithm in n^2
-  // using bounding boxes and the Octree would certainly be better
-  // for huge models...
+  if(_factory) return _factory->addCylinder(this, p1, p2, radius);
+  return 0;
+}
+
+GEntity *GModel::addTorus(std::vector<double> p1, std::vector<double> p2, 
+                          double radius1, double radius2)
+{
+  if(_factory) return _factory->addTorus(this, p1, p2, radius1, radius2);
+  return 0;
+}
 
+GEntity *GModel::addBlock(std::vector<double> p1, std::vector<double> p2)
+{
+  if(_factory) return _factory->addBlock(this, p1, p2);
+  return 0;
+}
+
+GEntity *GModel::addCone(std::vector<double> p1, std::vector<double> p2,
+                         double radius1, double radius2)
+{
+  if(_factory) return _factory->addCone(this, p1, p2,radius1, radius2);
+  return 0;
+}
+
+GModel *GModel::computeBooleanUnion(GModel *tool, int createNewModel)
+{
+  if(_factory)
+    return _factory->computeBooleanUnion(this, tool, createNewModel);
+  return 0;
+}
+
+GModel *GModel::computeBooleanIntersection(GModel *tool, int createNewModel)
+{
+  if(_factory)
+    return _factory->computeBooleanIntersection(this, tool, createNewModel);
+  return 0;
+}
+
+GModel *GModel::computeBooleanDifference(GModel *tool, int createNewModel)
+{
+  if(_factory)
+    return _factory->computeBooleanDifference(this, tool, createNewModel);
+  return 0;
+}
+
+static void computeDuplicates(GModel *model, 
+                              std::multimap<GVertex*, GVertex*> &Unique2Duplicates,
+                              std::map<GVertex*, GVertex*> &Duplicates2Unique,
+                              const double &eps)
+{
+  // currently we use a greedy algorithm in n^2 (using bounding boxes
+  // and the Octree would certainly be better for huge models...)
   std::list<GVertex*> v;
-  
-  v.insert(v.begin(),model->firstVertex(),model->lastVertex());
+  v.insert(v.begin(), model->firstVertex(), model->lastVertex());
 
   while(!v.empty()){
     GVertex *pv = *v.begin();
@@ -1422,91 +1548,79 @@ static void ComputeDuplicates(GModel * model,
     for ( std::multimap<GVertex*,GVertex*>::iterator it = Unique2Duplicates.begin(); 
           it != Unique2Duplicates.end(); ++it ){
       GVertex *unique = it->first;
-      const double d = sqrt ((unique->x()-pv->x())*(unique->x()-pv->x())+
-			     (unique->y()-pv->y())*(unique->y()-pv->y())+
-			     (unique->z()-pv->z())*(unique->z()-pv->z()));
+      const double d = sqrt((unique->x() - pv->x()) * (unique->x() - pv->x()) +
+                            (unique->y() - pv->y()) * (unique->y() - pv->y()) +
+                            (unique->z() - pv->z()) * (unique->z() - pv->z()));
       if (d <= eps) {
 	found = true;
-	Unique2Duplicates.insert(std::make_pair(unique,pv));
+	Unique2Duplicates.insert(std::make_pair(unique, pv));
 	Duplicates2Unique[pv] = unique;
 	break;
       }
     }
     if (!found) {
-      Unique2Duplicates.insert(std::make_pair(pv,pv));
+      Unique2Duplicates.insert(std::make_pair(pv, pv));
       Duplicates2Unique[pv] = pv;
     }
   }  
 }
 
-static void glueVerticesInEdges(GModel * model,
-				std::multimap<GVertex*,GVertex*>  & Unique2Duplicates,
-				std::map<GVertex*,GVertex*> & Duplicates2Unique){
+static void glueVerticesInEdges(GModel *model, 
+                                std::multimap<GVertex*, GVertex*> &Unique2Duplicates,
+				std::map<GVertex*, GVertex*> &Duplicates2Unique)
+{
   Msg::Debug("Gluing Edges");
-  for (GModel::eiter it = model->firstEdge(); it != model->lastEdge();++it){
+  for (GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
     GEdge *e = *it;
     GVertex *v1 = e->getBeginVertex();
     GVertex *v2 = e->getEndVertex();
     GVertex *replacementOfv1 = Duplicates2Unique[v1];
     GVertex *replacementOfv2 = Duplicates2Unique[v2];
     if ((v1 != replacementOfv1) || (v2 != replacementOfv2)){
-      Msg::Debug("Model Edge %d is re-build",e->tag());
-      e->replaceEndingPoints (replacementOfv1,replacementOfv2);
+      Msg::Debug("Model Edge %d is re-build", e->tag());
+      e->replaceEndingPoints (replacementOfv1, replacementOfv2);
     }
   }
 }
 
-static void ComputeDuplicates(GModel * model, 
-                              std::multimap<GEdge*,GEdge*>  & Unique2Duplicates,
-                              std::map<GEdge*,GEdge*> & Duplicates2Unique,
+static void computeDuplicates(GModel *model,
+                              std::multimap<GEdge*, GEdge*> &Unique2Duplicates,
+                              std::map<GEdge*,GEdge*> &Duplicates2Unique,
                               const double &eps)
 {
-  // in a first time, we use a greedy algorithm in n^2
-
-  // first check edges that have same endpoints
-
   std::list<GEdge*> e;
-  
-  e.insert(e.begin(),model->firstEdge(),model->lastEdge());
+  e.insert(e.begin(), model->firstEdge(), model->lastEdge());
 
   while(!e.empty()){
     GEdge *pe = *e.begin();
     // compute a point
     Range<double> r = pe->parBounds(0);
-    GPoint gp = pe->point(0.5*(r.low()+r.high()));
+    GPoint gp = pe->point(0.5 * (r.low() + r.high()));
     e.erase(e.begin());
     bool found = false;
     for (std::multimap<GEdge*,GEdge*>::iterator it = Unique2Duplicates.begin(); 
          it != Unique2Duplicates.end(); ++it ){
       GEdge *unique = it->first;
-            
-      //      printf ("checking %d %d\n",unique->tag(),pe->tag());
       // first check edges that have same endpoints
       if ((unique->getBeginVertex() == pe->getBeginVertex() && 
 	   unique->getEndVertex() == pe->getEndVertex()) ||
 	  (unique->getEndVertex() == pe->getBeginVertex() && 
 	   unique->getBeginVertex() == pe->getEndVertex())){
-	//	printf ("topology ok\n");
-	if (unique->geomType() == GEntity::Line &&
-	    pe->geomType() == GEntity::Line ){
+	if (unique->geomType() == GEntity::Line && pe->geomType() == GEntity::Line){
 	  found = true;
 	  Unique2Duplicates.insert(std::make_pair(unique,pe));
 	  Duplicates2Unique[pe] = unique;
-	  //	  printf ("D2U[%d] = %d\n",pe->tag(),unique->tag());
 	  break;	  
 	} 
 	double t;
 	GPoint gp2 = pe->closestPoint(SPoint3(gp.x(),gp.y(),gp.z()),t);
-	const double d = sqrt((gp.x()-gp2.x())*(gp.x()-gp2.x())+
-			      (gp.y()-gp2.y())*(gp.y()-gp2.y())+
-			      (gp.z()-gp2.z())*(gp.z()-gp2.z()));
-	//	printf ("closest point @ distance %g\n",d);
+	const double d = sqrt((gp.x() - gp2.x()) * (gp.x() - gp2.x()) +
+			      (gp.y() - gp2.y()) * (gp.y() - gp2.y()) +
+			      (gp.z() - gp2.z()) * (gp.z() - gp2.z()));
 	if (t >= r.low() && t <= r.high() && d <= eps) {
-	  //	  printf ("both are close -> glue\n");
 	  found = true;
 	  Unique2Duplicates.insert(std::make_pair(unique,pe));
 	  Duplicates2Unique[pe] = unique;
-	  //	  printf ("D2U[%d] = %d\n",pe->tag(),unique->tag());
 	  break;
 	}
       }
@@ -1514,48 +1628,39 @@ static void ComputeDuplicates(GModel * model,
     if (!found) {
       Unique2Duplicates.insert(std::make_pair(pe,pe));
       Duplicates2Unique[pe] = pe;
-      //      printf ("D2U[%d] = %d\n",pe->tag(),pe->tag());
     }
   }
-
-  //  for(std::map<GEdge*,GEdge*>::iterator it = Duplicates2Unique.begin(); it != Duplicates2Unique.end(); ++it){
-  //    printf("E(%d) -- E(%d)\n",it->first->tag(),it->second->tag());
-  //  }
-
 }
 
-static void glueEdgesInFaces(GModel * model,
-			     std::multimap<GEdge*,GEdge*>  & Unique2Duplicates,
-			     std::map<GEdge*,GEdge*> & Duplicates2Unique)
+static void glueEdgesInFaces(GModel *model, 
+                             std::multimap<GEdge*, GEdge*> &Unique2Duplicates,
+			     std::map<GEdge*, GEdge*> &Duplicates2Unique)
 {
   Msg::Debug("Gluing Model Faces");
-  for (GModel::fiter it = model->firstFace(); it != model->lastFace();++it){
+  for (GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){
     GFace *f = *it;
     bool aDifferenceExists = false;
     std::list<GEdge*> old = f->edges(), enew;
-    //    printf("face %d %d edges\n",(*it)->tag(),old.size());
-    for (std::list<GEdge*>::iterator eit = old.begin(); eit !=old.end() ; eit++){
+    for (std::list<GEdge*>::iterator eit = old.begin(); eit !=old.end(); eit++){
       GEdge *temp = Duplicates2Unique[*eit];
       enew.push_back(temp);
-      //      printf("edge %d vs. %d\n",(*eit)->tag(),temp->tag());
       if (temp != *eit){
 	aDifferenceExists = true;
       }
     }
     if (aDifferenceExists){
-      Msg::Debug("Model Face %d is re-build",f->tag());
+      Msg::Debug("Model Face %d is re-build", f->tag());
       f->replaceEdges (enew);
     }
   }
 }
 
-static void ComputeDuplicates(GModel * model, 
-                              std::multimap<GFace*,GFace*>  & Unique2Duplicates,
-                              std::map<GFace*,GFace*> & Duplicates2Unique,
+static void computeDuplicates(GModel *model, 
+                              std::multimap<GFace*, GFace*> &Unique2Duplicates,
+                              std::map<GFace*,GFace*> &Duplicates2Unique, 
                               const double &eps)
 {
   std::list<GFace*> f;
-  
   f.insert(f.begin(),model->firstFace(),model->lastFace());
 
   while(!f.empty()){
@@ -1564,12 +1669,13 @@ static void ComputeDuplicates(GModel * model,
     Range<double> r = pf->parBounds(0);
     Range<double> s = pf->parBounds(1);
     // FIXME : this is WRONG (the point can be in a hole ;-) 
-    GPoint gp = pf->point(SPoint2(0.5*(r.low()+r.high()),0.5*(s.low()+s.high())));
+    GPoint gp = pf->point(SPoint2(0.5 * (r.low() + r.high()), 0.5 * (s.low() + s.high())));
     f.erase(f.begin());
     std::list<GEdge*> pf_edges = pf->edges();
     pf_edges.sort();
     bool found = false;
-    for ( std::multimap<GFace*,GFace*>::iterator it = Unique2Duplicates.begin(); it != Unique2Duplicates.end() ; ++it ){
+    for (std::multimap<GFace*,GFace*>::iterator it = Unique2Duplicates.begin();
+         it != Unique2Duplicates.end(); ++it){
       GFace *unique = it->first;      
       std::list<GEdge*> unique_edges = unique->edges();
       if (unique_edges.size() == pf_edges.size()){
@@ -1582,16 +1688,13 @@ static void ComputeDuplicates(GModel * model,
 	  if (*it_pf != *it_unique) all_similar = false;
 	}
 	if (all_similar){
-	  if (unique->geomType() == GEntity::Plane &&
-	      pf->geomType() == GEntity::Plane ){
+	  if (unique->geomType() == GEntity::Plane && pf->geomType() == GEntity::Plane){
 	    found = true;
 	    Unique2Duplicates.insert(std::make_pair(unique,pf));
 	    Duplicates2Unique[pf] = unique;
-	    //	    printf ("F: D2U[%d] = %d\n",pf->tag(),unique->tag());
 	    break;	  
 	  } 
 	  double t[2]={0,0};
-	  //	  GPoint gp2 = pf->closestPoint(SPoint3(gp.x(),gp.y(),gp.z()),t);
 	  const double d = 1.0;
 	  /*sqrt((gp.x()-gp2.x())*(gp.x()-gp2.x())+
 	    (gp.y()-gp2.y())*(gp.y()-gp2.y())+
@@ -1599,11 +1702,9 @@ static void ComputeDuplicates(GModel * model,
 	  //	printf ("closest point @ distance %g\n",d);
 	  if (t[0] >= r.low() && t[0] <= r.high() && 
 	      t[1] >= s.low() && t[1] <= s.high() && d <= eps) {
-	    //	  printf ("both are close -> glue\n");
 	    found = true;
 	    Unique2Duplicates.insert(std::make_pair(unique,pf));
 	    Duplicates2Unique[pf] = unique;
-	    //	  printf ("F : D2U[%d] = %d\n",pf->tag(),unique->tag());
 	    break;
 	  }
 	}
@@ -1612,34 +1713,28 @@ static void ComputeDuplicates(GModel * model,
     if (!found) {
       Unique2Duplicates.insert(std::make_pair(pf,pf));
       Duplicates2Unique[pf] = pf;
-      //      printf ("F : D2U[%d] = %d\n",pf->tag(),pf->tag());
     }
   }
-  
-  //  for(std::map<GFace*,GFace*>::iterator it = Duplicates2Unique.begin(); it != Duplicates2Unique.end(); ++it){
-  //    printf("F(%d) -- F(%d)\n",it->first->tag(),it->second->tag());
-  //  }  
 }
 
-static void glueFacesInRegions(GModel * model,
-			       std::multimap<GFace*,GFace*>  & Unique2Duplicates,
-			       std::map<GFace*,GFace*> & Duplicates2Unique)
+static void glueFacesInRegions(GModel *model,
+                               std::multimap<GFace*, GFace*> &Unique2Duplicates,
+			       std::map<GFace*, GFace*> &Duplicates2Unique)
 {
   Msg::Debug("Gluing Regions");
   for (GModel::riter it = model->firstRegion(); it != model->lastRegion();++it){
     GRegion *r = *it;
     bool aDifferenceExists = false;
     std::list<GFace*> old = r->faces(), fnew;
-    for (std::list<GFace*>::iterator fit = old.begin(); fit !=old.end() ; fit++){
+    for (std::list<GFace*>::iterator fit = old.begin(); fit != old.end(); fit++){
       GFace *temp = Duplicates2Unique[*fit];
       fnew.push_back(temp);
-      //      printf("edge %d vs. %d\n",(*eit)->tag(),temp->tag());
       if (temp != *fit){
 	aDifferenceExists = true;
       }
     }
     if (aDifferenceExists){
-      Msg::Debug("Model Region %d is re-build",r->tag());
+      Msg::Debug("Model Region %d is re-build", r->tag());
       r->replaceFaces (fnew);
     }
   }
@@ -1649,21 +1744,21 @@ void GModel::glue(const double &eps)
 {
   {
     std::multimap<GVertex*,GVertex*> Unique2Duplicates;
-    std::map<GVertex*,GVertex*>      Duplicates2Unique;
-    ComputeDuplicates  (this, Unique2Duplicates,Duplicates2Unique, eps);
-    glueVerticesInEdges(this, Unique2Duplicates,Duplicates2Unique);
+    std::map<GVertex*,GVertex*> Duplicates2Unique;
+    computeDuplicates(this, Unique2Duplicates, Duplicates2Unique, eps);
+    glueVerticesInEdges(this, Unique2Duplicates, Duplicates2Unique);
   }
   {
     std::multimap<GEdge*,GEdge*> Unique2Duplicates;
-    std::map<GEdge*,GEdge*>      Duplicates2Unique;
-    ComputeDuplicates  (this, Unique2Duplicates,Duplicates2Unique, eps);
-    glueEdgesInFaces(this, Unique2Duplicates,Duplicates2Unique);
+    std::map<GEdge*,GEdge*> Duplicates2Unique;
+    computeDuplicates(this, Unique2Duplicates, Duplicates2Unique, eps);
+    glueEdgesInFaces(this, Unique2Duplicates, Duplicates2Unique);
   }    
   {
     std::multimap<GFace*,GFace*> Unique2Duplicates;
-    std::map<GFace*,GFace*>      Duplicates2Unique;
-    ComputeDuplicates  (this, Unique2Duplicates,Duplicates2Unique, eps);
-    glueFacesInRegions(this, Unique2Duplicates,Duplicates2Unique);
+    std::map<GFace*,GFace*> Duplicates2Unique;
+    computeDuplicates(this, Unique2Duplicates, Duplicates2Unique, eps);
+    glueFacesInRegions(this, Unique2Duplicates, Duplicates2Unique);
   }    
 }
 
@@ -1673,6 +1768,7 @@ void GModel::registerBindings(binding *b)
 {
   classBinding *cb = b->addClass<GModel>("GModel");
   cb->setDescription("A GModel contains a geometry and its mesh.");
+
   methodBinding *cm;
   cm = cb->addMethod("mesh", &GModel::mesh);
   cm->setArgNames("dim", NULL);
@@ -1711,13 +1807,72 @@ void GModel::registerBindings(binding *b)
   cm->setArgNames("tag", NULL);
   cm = cb->addMethod("getVertexByTag", &GModel::getVertexByTag);
   cm->setDescription("access a geometrical vertex by tag");
-  cm->setArgNames("tag",NULL);
+  cm->setArgNames("tag", NULL);
   cm = cb->addMethod("getRegionByTag", &GModel::getRegionByTag);
   cm->setDescription("access a geometrical region by tag");
   cm->setArgNames("tag", NULL);
+
+  cm = cb->addMethod("addVertex", &GModel::addVertex);
+  cm->setDescription("create a new vertex at position (x, y, z), with target "
+                     "mesh size lc");
+  cm->setArgNames("x", "y", "z", "lc", NULL);
+  cm = cb->addMethod("addLine", &GModel::addLine);
+  cm->setDescription("create a straight line going from v1 to v2");
+  cm->setArgNames("v1", "v2", NULL);
+  cm = cb->addMethod("addBezier", &GModel::addBezier);
+  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("addCircleArcCenter", &GModel::addCircleArcCenter);
+  cm->setDescription("create a circle arc going from v1 to v2 with its center "
+                     "at (x,y,z)");
+  cm->setArgNames("x", "y", "z", "v1", "v2", NULL);
+  cm = cb->addMethod("addCircleArc3Points", &GModel::addCircleArc3Points);
+  cm->setDescription("create a circle arc going from v1 to v2 with an "
+                     "intermediary point Xi(x,y,z)");
+  cm->setArgNames("x", "y", "z", "v1", "v2", NULL);
+  cm = cb->addMethod("revolve", &GModel::revolve);
+  cm->setDescription("revolve an entity of a given angle. Axis is defined by 2 "
+                     "points "
+                     "in a full Matrix(2,3)");
+  cm->setArgNames("entity", "angle", "axis", NULL);
+  cm = cb->addMethod("extrude", &GModel::extrude);
+  cm->setDescription("extrudes an entity. Axis is defined by 2 points in a full "
+                     "Matrix(2,3)");
+  cm->setArgNames("entity", "axis", NULL);
+  cm = cb->addMethod("addSphere", &GModel::addSphere);
+  cm->setDescription("add a sphere");
+  cm->setArgNames("xc", "yc", "zc", "radius", NULL);
+  cm = cb->addMethod("addBlock", &GModel::addBlock);
+  cm->setDescription("add a block");
+  cm->setArgNames("{x1,y1,z1}", "{x2,y2,z2}", NULL);
+  cm = cb->addMethod("addCone", &GModel::addCone);
+  cm->setDescription("add a cone");
+  cm->setArgNames("{x1,y1,z1}","{x2,y2,z2}","R1","R2",NULL);
+  cm = cb->addMethod("addCylinder", &GModel::addCylinder);
+  cm->setDescription("add a cylinder");
+  cm->setArgNames("{x1,y1,z1}","{x2,y2,z2}", "R",NULL);
+  cm = cb->addMethod("addTorus", &GModel::addTorus);
+  cm->setDescription("add a torus");
+  cm->setArgNames("{x1,y1,z1}","{x2,y2,z2}","R1","R2",NULL);
+  cm = cb->addMethod("computeUnion", &GModel::computeBooleanUnion);
+  cm->setDescription("compute the boolean union of the model with another one "
+                     "(tool). The third parameter tells if a new model has to "
+                     "be created");
+  cm->setArgNames("tool", "createNewGModel",NULL);
+  cm = cb->addMethod("computeIntersection", &GModel::computeBooleanIntersection);
+  cm->setDescription("compute the boolean intersection of the model with another "
+                     "one. The third parameter tells if a new model has to be created");
+  cm->setArgNames("tool","createNewGModel",NULL);
+  cm = cb->addMethod("computeDifference", &GModel::computeBooleanDifference);
+  cm->setDescription("compute the boolean difference of the model with another "
+                     "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->setConstructor<GModel>();
   cm->setDescription("Create an empty GModel");
 }
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 93b0b865c5..9d70d4563c 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -17,7 +17,7 @@
 #include "GRegion.h"
 #include "SPoint3.h"
 #include "SBoundingBox3d.h"
-#include "discreteFace.h"
+#include "fullMatrix.h"
 
 class Octree;
 class FM_Internals;
@@ -29,7 +29,7 @@ class CGNSOptions;
 class gLevelset;
 class discreteFace;
 class binding;
-class OCCFactory;
+class GModelFactory;
 
 // A geometric model. The model is a "not yet" non-manifold B-Rep.
 class GModel
@@ -73,6 +73,9 @@ class GModel
   void _createFMInternals();
   void _deleteFMInternals();
 
+  // CAD creation factory
+  GModelFactory *_factory;
+
   // characteristic length (mesh size) fields
   FieldManager *_fields;
 
@@ -342,13 +345,38 @@ 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 (const double &eps);
+  // change the entity creation factory
+  void setFactory(std::string name);
+
+  // create brep geometry entities using the factory
+  GVertex *addVertex(double x, double y, double z, double lc);
+  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);
+  GEntity *revolve(GEntity *e, double angle, fullMatrix<double> *axis);
+  GEntity *extrude(GEntity *e, fullMatrix<double> *axis);
+
+  // create solid geometry primitives using the factory
+  GEntity *addSphere(double cx, double cy, double cz, double radius);
+  GEntity *addCylinder(std::vector<double> p1, std::vector<double> p2, double radius);
+  GEntity *addTorus(std::vector<double> p1, std::vector<double> p2, double radius1,
+                    double radius2);
+  GEntity *addBlock(std::vector<double> p1, std::vector<double> p2);
+  GEntity *addCone(std::vector<double> p1, std::vector<double> p2, double radius1,
+                   double radius2);
+
+  // boolean operators acting on 2 models
+  GModel *computeBooleanUnion(GModel *tool, int createNewModel);
+  GModel *computeBooleanIntersection(GModel *tool, int createNewModel);
+  GModel *computeBooleanDifference(GModel *tool, int createNewModel);
+
+  // glue entities in the model (i.e., assume a tolerance eps and
+  // merge vertices that are too close, then merge edges, faces and
+  // regions). Warning: 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
+  void glue(const double &eps);
 
   // build a new GModel by cutting the elements crossed by the levelset ls
   // if cutElem is set to false, split the model without cutting the elements
@@ -365,6 +393,15 @@ class GModel
                               std::vector<int> &elementary,
                               std::vector<int> &partition);
 
+  // Store mesh elements of a chain in a new elementary and physical entity
+  void storeChain(int dim, std::map<int, std::vector<MElement*> > &entityMap,
+                  std::map<int, std::map<int, std::string> > &physicalMap)
+  {
+    _storeElementsInEntities(entityMap);
+    _storePhysicalTagsInEntities(dim, physicalMap);
+    _associateEntityWithMeshVertices();
+  }
+
   // Gmsh native CAD format (readGEO is static, since it can create
   // multiple models)
   static int readGEO(const std::string &name);
@@ -460,20 +497,10 @@ class GModel
   int writeDIFF(const std::string &name, bool binary=false,
                bool saveAll=false, double scalingFactor=1.0);
 
-
   void save(std::string fileName);
   void load(std::string fileName);
 
   static void registerBindings(binding *b);
-
-  // Store mesh elements of a chain in a new elementary and physical entity
-  void storeChain(int dim, std::map<int, std::vector<MElement*> > &entityMap,
-                  std::map<int, std::map<int, std::string> > &physicalMap)
-  {
-    _storeElementsInEntities(entityMap);
-    _storePhysicalTagsInEntities(dim, physicalMap);
-    _associateEntityWithMeshVertices();
-  }
 };
 
 #endif
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 649d51285f..1ef30f0b8d 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -1,6 +1,12 @@
+// 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 "GModelFactory.h"
 
 #if defined(HAVE_OCC)
+
 #include "OCCIncludes.h"
 #include "GModelIO_OCC.h"
 #include "OCCVertex.h"
@@ -20,27 +26,23 @@
 #include <Geom_TrimmedCurve.hxx>
 #include <Geom_BezierCurve.hxx>
 
-OCCFactory::OCCFactory() : GModelFactory()
-{
-}
-
-GVertex * OCCFactory::createVertex (GModel *_gm, double x, double y, double z, double lc)
+GVertex *OCCFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
 
   gp_Pnt aPnt;
   aPnt = gp_Pnt(x, y, z);
-  BRepBuilderAPI_MakeVertex mkVertex (aPnt);
-  TopoDS_Vertex occv =  mkVertex.Vertex();
+  BRepBuilderAPI_MakeVertex mkVertex(aPnt);
+  TopoDS_Vertex occv = mkVertex.Vertex();
   
-  return _gm->_occ_internals->addVertexToModel(_gm,occv);
+  return gm->_occ_internals->addVertexToModel(gm, occv);
 }
 
-GEdge * OCCFactory::createLine (GModel *_gm, GVertex *start, GVertex *end)
+GEdge *OCCFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
 
   OCCVertex *occv1 = dynamic_cast<OCCVertex*>(start);
   OCCVertex *occv2 = dynamic_cast<OCCVertex*>(end);
@@ -54,249 +56,256 @@ GEdge * OCCFactory::createLine (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, start, end);
 }
 
-GEdge *OCCFactory::createCircleArc (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;
+  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());
+  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());
 
   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);
+    return gm->_occ_internals->addEdgeToModel(gm, occEdge, start, end);
   }
   else if (method == GModelFactory::CENTER_START_END){
     Standard_Real Radius = aP1.Distance(aP2);
-    gce_MakeCirc MC(aP1,gce_MakePln(aP1, aP2, aP3).Value(),Radius);
+    gce_MakeCirc MC(aP1,gce_MakePln(aP1, aP2, aP3).Value(), Radius);
     const gp_Circ& Circ = MC.Value();
-    Standard_Real Alpha1 = ElCLib::Parameter(Circ,aP1);
-    Standard_Real Alpha2 = ElCLib::Parameter(Circ,aP3);
+    Standard_Real Alpha1 = ElCLib::Parameter(Circ, aP1);
+    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);
+    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);
+    return gm->_occ_internals->addEdgeToModel(gm, occEdge, start, end);
   }
   return 0;
 }
 
-GEdge *OCCFactory::createSpline (GModel *_gm, const splineType &type,
-				 GVertex *start, GVertex *end, 
-				 fullMatrix<double> *points)
+GEdge *OCCFactory::addSpline(GModel *gm, const splineType &type,
+                             GVertex *start, GVertex *end, 
+                             fullMatrix<double> *points)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
 
   OCCEdge *occEd = 0;
   int nbControlPoints = points->size1();
-  TColgp_Array1OfPnt ctrlPoints (1, nbControlPoints+2);
+  TColgp_Array1OfPnt ctrlPoints (1, nbControlPoints + 2);
   int index = 1;
-  ctrlPoints.SetValue(index++,gp_Pnt(start->x(),start->y(),start->z()));  
+  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()));  
+  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);
+    return gm->_occ_internals->addEdgeToModel(gm, bez, start, end);
   } 
-  else {
-    throw;
-  }
+  Msg::Error("Non-bezier splines not implemented yet");
   return 0;
 }
 
-GEntity *OCCFactory::revolve (GModel *_gm, GEntity* base, 
-			      double x1, double y1, double z1, 
-			      double x2, double y2, double z2,
-			      double angle)
+GEntity *OCCFactory::revolve(GModel *gm, GEntity* base, 
+                             double x1, double y1, double z1, 
+                             double x2, double y2, double z2,
+                             double angle)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
 
-  gp_Dir direction (x2-x1,y2-y1,z2-z1);
-  gp_Ax1 axisOfRevolution (gp_Pnt (x1,y1,z1),direction);
-  BRepPrimAPI_MakeRevol MR (*(TopoDS_Shape*)base->getNativePtr(), 
-			    axisOfRevolution, angle, Standard_False);
+  gp_Dir direction(x2 - x1, y2 - y1, z2 - z1);
+  gp_Ax1 axisOfRevolution(gp_Pnt(x1, y1, z1), direction);
+  BRepPrimAPI_MakeRevol MR(*(TopoDS_Shape*)base->getNativePtr(), 
+                           axisOfRevolution, angle, Standard_False);
   GEntity *ret;
   if (base->cast2Vertex()){
-    TopoDS_Edge result =  TopoDS::Edge(MR.Shape());
-    ret  = _gm->_occ_internals->addEdgeToModel(_gm,result);
+    TopoDS_Edge result = TopoDS::Edge(MR.Shape());
+    ret = gm->_occ_internals->addEdgeToModel(gm, result);
   }
   if (base->cast2Edge()){
-    TopoDS_Face result =  TopoDS::Face(MR.Shape());
-    ret =  _gm->_occ_internals->addFaceToModel(_gm,result);    
+    TopoDS_Face result = TopoDS::Face(MR.Shape());
+    ret = gm->_occ_internals->addFaceToModel(gm, result);
   }
   if (base->cast2Face()){
-    TopoDS_Solid result =  TopoDS::Solid(MR.Shape());
-    ret =  _gm->_occ_internals->addRegionToModel(_gm,result);    
+    TopoDS_Solid result = TopoDS::Solid(MR.Shape());
+    ret = gm->_occ_internals->addRegionToModel(gm, result);
   }
-  _gm->glue(Precision::Confusion());
+  gm->glue(Precision::Confusion());
   return ret;
 }
 
-GEntity *OCCFactory::extrude(GModel *_gm, GEntity* base, 
+GEntity *OCCFactory::extrude(GModel *gm, GEntity* base, 
                              double x1, double y1, double z1, 
                              double x2, double y2, double z2)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
 
-  gp_Vec direction (gp_Pnt (x1,y1,z1),gp_Pnt (x2,y2,z2));
-  gp_Ax1 axisOfRevolution (gp_Pnt (x1,y1,z1),direction);
+  gp_Vec direction(gp_Pnt(x1, y1, z1), gp_Pnt(x2, y2, z2));
+  gp_Ax1 axisOfRevolution(gp_Pnt(x1, y1, z1), direction);
 
-  BRepPrimAPI_MakePrism MP(*(TopoDS_Shape*)base->getNativePtr(), direction, Standard_False);
+  BRepPrimAPI_MakePrism MP(*(TopoDS_Shape*)base->getNativePtr(), direction, 
+                           Standard_False);
 
   GEntity *ret;
   
   if (base->cast2Vertex()){
-    TopoDS_Edge result =  TopoDS::Edge(MP.Shape());
-    ret =  _gm->_occ_internals->addEdgeToModel(_gm,result);
+    TopoDS_Edge result = TopoDS::Edge(MP.Shape());
+    ret = gm->_occ_internals->addEdgeToModel(gm, result);
   }
   if (base->cast2Edge()){
-    TopoDS_Face result =  TopoDS::Face(MP.Shape());
-    ret = _gm->_occ_internals->addFaceToModel(_gm,result);    
+    TopoDS_Face result = TopoDS::Face(MP.Shape());
+    ret = gm->_occ_internals->addFaceToModel(gm, result);    
   }
   if (base->cast2Face()){
-    TopoDS_Solid result =  TopoDS::Solid(MP.Shape());
-    ret = _gm->_occ_internals->addRegionToModel(_gm,result);    
+    TopoDS_Solid result = TopoDS::Solid(MP.Shape());
+    ret = gm->_occ_internals->addRegionToModel(gm, result);    
   }
-  _gm->glue(Precision::Confusion());
+  gm->glue(Precision::Confusion());
   return ret;
 }
 
-GEntity *OCCFactory::sphere (GModel *_gm,  double xc, double yc, double zc, double radius)
+GEntity *OCCFactory::addSphere(GModel *gm, double xc, double yc, double zc, double radius)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
-
-  gp_Pnt aP(xc,yc,zc);  
-  _gm->_occ_internals->buildShapeFromLists(BRepPrimAPI_MakeSphere(aP, radius).Shape());
-  _gm->destroy();
-  _gm->_occ_internals->buildLists();
-  _gm->_occ_internals->buildGModel(_gm);  
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
+
+  gp_Pnt aP(xc, yc, zc);  
+  gm->_occ_internals->buildShapeFromLists(BRepPrimAPI_MakeSphere(aP, radius).Shape());
+  gm->destroy();
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);  
+  return 0; // FIXME: should return GRegion?
 }
 
-GEntity *OCCFactory::cylinder(GModel *_gm,  std::vector<double> p1, std::vector<double> p2,
-                              double radius)
+GEntity *OCCFactory::addCylinder(GModel *gm, std::vector<double> p1,
+                                 std::vector<double> p2, double radius)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
-
-  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]; 
-
-  gp_Pnt aP(x1,y1,z1);
-  const double H = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)); 
-  gp_Vec aV((x2-x1)/H,(y2-y1)/H,(z2-z1)/H);
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
+
+  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];
+
+  gp_Pnt aP(x1, y1, z1);
+  const double H = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + 
+                        (z2 - z1) * (z2 - z1));
+  gp_Vec aV((x2 - x1) / H, (y2 - y1) / H, (z2 - z1) / H);
   gp_Ax2 anAxes(aP, aV);
-  BRepPrimAPI_MakeCylinder MC (anAxes, radius, H);
+  BRepPrimAPI_MakeCylinder MC(anAxes, radius, H);
   MC.Build();
   if (!MC.IsDone()) {
     Msg::Error("Cylinder can't be computed from the given parameters");
     return 0;
   }
-  _gm->_occ_internals->buildShapeFromLists(MC.Shape());
-  _gm->destroy();
-  _gm->_occ_internals->buildLists();
-  _gm->_occ_internals->buildGModel(_gm);  
-  return 0;
+  gm->_occ_internals->buildShapeFromLists(MC.Shape());
+  gm->destroy();
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);
+  return 0;  // FIXME: should return GRegion?
 }
 
-GEntity *OCCFactory::torus(GModel *_gm,  std::vector<double> p1, std::vector<double> p2,
-                           double radius1, double radius2)
+GEntity *OCCFactory::addTorus(GModel *gm, std::vector<double> p1, 
+                              std::vector<double> p2, double radius1, 
+                              double radius2)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
-
-  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]; 
-  gp_Pnt aP(x1,y1,z1);
-  const double H = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)); 
-  gp_Vec aV((x2-x1)/H,(y2-y1)/H,(z2-z1)/H);
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
+
+  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];
+  gp_Pnt aP(x1, y1, z1);
+  const double H = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) +
+                        (z2 - z1) * (z2 - z1)); 
+  gp_Vec aV((x2 - x1) / H, (y2 - y1) / H, (z2 - z1) / H);
   gp_Ax2 anAxes(aP, aV);
-  BRepPrimAPI_MakeTorus MC (anAxes, radius1, radius2);
+  BRepPrimAPI_MakeTorus MC(anAxes, radius1, radius2);
   MC.Build();
   if (!MC.IsDone()) {
     Msg::Error("Cylinder can't be computed from the given parameters");
     return 0;
   }
-  _gm->_occ_internals->buildShapeFromLists(MC.Shape());
-  _gm->destroy();
-  _gm->_occ_internals->buildLists();
-  _gm->_occ_internals->buildGModel(_gm);  
-  return 0;
+  gm->_occ_internals->buildShapeFromLists(MC.Shape());
+  gm->destroy();
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);
+  return 0;  // FIXME: should return GRegion?
 }
 
-GEntity *OCCFactory::cone (GModel *_gm,  std::vector<double> p1, std::vector<double> p2, 
-                           double radius1, double radius2)
+GEntity *OCCFactory::addCone(GModel *gm,  std::vector<double> p1, 
+                             std::vector<double> p2, double radius1, 
+                             double radius2)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
-
-  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]; 
-
-  gp_Pnt aP(x1,y1,z1);
-  const double H = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)); 
-  gp_Vec aV((x2-x1)/H,(y2-y1)/H,(z2-z1)/H);
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
+
+  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];
+
+  gp_Pnt aP(x1, y1, z1);
+  const double H = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) +
+                        (z2 - z1) * (z2 - z1)); 
+  gp_Vec aV((x2 - x1) / H, (y2 - y1) / H, (z2 - z1) / H);
   gp_Ax2 anAxes(aP, aV);
-  BRepPrimAPI_MakeCone MC (anAxes, radius1, radius2, H);
+  BRepPrimAPI_MakeCone MC(anAxes, radius1, radius2, H);
   MC.Build();
   if (!MC.IsDone()) {
     Msg::Error("Cylinder can't be computed from the given parameters");
     return 0;
   }
-  _gm->_occ_internals->buildShapeFromLists(MC.Shape());
-  _gm->destroy();
-  _gm->_occ_internals->buildLists();
-  _gm->_occ_internals->buildGModel(_gm);  
-  return 0;
+  gm->_occ_internals->buildShapeFromLists(MC.Shape());
+  gm->destroy();
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);  
+  return 0; // FIXME: should return GRegion?
 }
 
-GEntity *OCCFactory::block  (GModel *_gm, std::vector<double> p1, std::vector<double> p2)
+GEntity *OCCFactory::addBlock(GModel *gm, std::vector<double> p1, 
+                              std::vector<double> p2)
 {
-  if (!_gm->_occ_internals)
-    _gm->_occ_internals = new OCC_Internals;
+  if (!gm->_occ_internals)
+    gm->_occ_internals = new OCC_Internals;
 
-  gp_Pnt P1(p1[0],p1[1],p1[2]);
-  gp_Pnt P2(p2[0],p2[1],p2[2]);
+  gp_Pnt P1(p1[0], p1[1], p1[2]);
+  gp_Pnt P2(p2[0], p2[1], p2[2]);
   BRepPrimAPI_MakeBox MB(P1, P2);
   MB.Build();
   if (!MB.IsDone()) {
     Msg::Error("Box can not be computed from the given point");
     return 0;
   }
-  _gm->_occ_internals->buildShapeFromLists(MB.Shape());
-  _gm->destroy();
-  _gm->_occ_internals->buildLists();
-  _gm->_occ_internals->buildGModel(_gm);  
+  gm->_occ_internals->buildShapeFromLists(MB.Shape());
+  gm->destroy();
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);
+  return 0;  // FIXME: should return GRegion?
 }
 
-GModel *OCCFactory::booleanUnion (GModel* obj, GModel* tool, int createNewModel)
+GModel *OCCFactory::computeBooleanUnion(GModel *obj, GModel *tool, int createNewModel)
 {
-  OCC_Internals *occ_obj  = obj->getOCCInternals();
+  OCC_Internals *occ_obj = obj->getOCCInternals();
   OCC_Internals *occ_tool = tool->getOCCInternals();
 
   if (!occ_obj || !occ_tool)return NULL;
@@ -307,17 +316,17 @@ GModel *OCCFactory::booleanUnion (GModel* obj, GModel* tool, int createNewModel)
     temp->_occ_internals->addShapeToLists(occ_obj->getShape());
     obj = temp;
   }
-  obj->_occ_internals->applyBooleanOperator(occ_tool->getShape(),OCC_Internals::Fuse);
+  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::booleanDifference (GModel* obj, GModel* tool, int createNewModel)
+GModel *OCCFactory::computeBooleanDifference(GModel *obj, GModel *tool,
+                                             int createNewModel)
 {
-  OCC_Internals *occ_obj  = obj->getOCCInternals();
+  OCC_Internals *occ_obj = obj->getOCCInternals();
   OCC_Internals *occ_tool = tool->getOCCInternals();
 
   if (!occ_obj || !occ_tool)return NULL;
@@ -328,19 +337,20 @@ GModel *OCCFactory::booleanDifference (GModel* obj, GModel* tool, int createNewM
     temp->_occ_internals->addShapeToLists(occ_obj->getShape());
     obj = temp;
   }
-  obj->getOCCInternals()->applyBooleanOperator(occ_tool->getShape(),OCC_Internals::Cut);
+  obj->getOCCInternals()->applyBooleanOperator(occ_tool->getShape(), 
+                                               OCC_Internals::Cut);
   obj->destroy();
   obj->_occ_internals->buildLists();
   obj->_occ_internals->buildGModel(obj);
   return obj;
 }
 
-GModel *OCCFactory::booleanIntersection (GModel* obj, GModel* tool, int createNewModel)
+GModel *OCCFactory::computeBooleanIntersection(GModel *obj, GModel* tool, 
+                                               int createNewModel)
 {
-  OCC_Internals *occ_obj  = obj->getOCCInternals();
+  OCC_Internals *occ_obj = obj->getOCCInternals();
   OCC_Internals *occ_tool = tool->getOCCInternals();
 
-
   if (!occ_obj || !occ_tool)return NULL;
 
   if (createNewModel){
@@ -358,71 +368,3 @@ GModel *OCCFactory::booleanIntersection (GModel* obj, GModel* tool, int createNe
 }
 
 #endif
-
-#include "Bindings.h"
-void GModelFactory::registerBindings (binding *b)
-{
-  classBinding *cb = b->addClass<GModelFactory>("GModelFactory");
-  cb->setDescription("an interface to Gentity Construction");
-  methodBinding *mb;
-  mb = cb->addMethod("createVertex", &GModelFactory::createVertex);
-  mb->setDescription("creates a GVertex");
-  mb->setArgNames("model", "x", "y","z","lc",NULL);
-  mb = cb->addMethod("createLine", &GModelFactory::createLine);
-  mb->setDescription("creates a Line going from v1 to v2");
-  mb->setArgNames("model","v1", "v2",NULL);
-  mb = cb->addMethod("createBezierLine", &GModelFactory::createBezier);
-  mb->setDescription("creates a Spline going from v1 to v2 and with some control points listed in a fullMatrix(N,3)");
-  mb->setArgNames("model","v1", "v2", "controlPoints", NULL);
-  mb = cb->addMethod("createCircleArcCenter", &GModelFactory::createCircleArcCenter);
-  mb->setDescription("creates a circle arc going from v1 to v2 with its center Xc(x,y,z)");
-  mb->setArgNames("model","x","y","z","v1", "v2", NULL);
-  mb = cb->addMethod("createCircleArc3Points", &GModelFactory::createCircleArc3Points);
-  mb->setDescription("creates a circle arc going from v1 to v2 with an intermediary point Xi(x,y,z)");
-  mb->setArgNames("model","x","y","z","v1", "v2", NULL);
-  mb = cb->addMethod("revolve", &GModelFactory::revolve_);
-  mb->setDescription("revolves an entity of an angle. Axis is defined by 2 points in a full Matrix(2,3)");
-  mb->setArgNames("model","entity", "angle", "axis",NULL);
-  mb = cb->addMethod("extrude", &GModelFactory::extrude_);
-  mb->setDescription("extrudes an entity. Axis is defined by 2 points in a full Matrix(2,3)");
-  mb->setArgNames("model","entity", "axis",NULL);
-
-  mb = cb->addMethod("sphere", &GModelFactory::sphere);
-  mb->setDescription("builds a sphere");
-  mb->setArgNames("model","xc","yc","zc","radius",NULL);
-
-  mb = cb->addMethod("block", &GModelFactory::block);
-  mb->setDescription("builds a block");
-  mb->setArgNames("model","{x1,y1,z1}","{x2,y2,z2}",NULL);
-
-  mb = cb->addMethod("cone", &GModelFactory::cone);
-  mb->setDescription("builds a cone");
-  mb->setArgNames("model","{x1,y1,z1}","{x2,y2,z2}","R1","R2",NULL);
-
-  mb = cb->addMethod("cylinder", &GModelFactory::cylinder);
-  mb->setDescription("builds a cylinder");
-  mb->setArgNames("model","{x1,y1,z1}","{x2,y2,z2}", "R",NULL);
-
-  mb = cb->addMethod("torus", &GModelFactory::torus);
-  mb->setDescription("builds a torus");
-  mb->setArgNames("model","{x1,y1,z1}","{x2,y2,z2}","R1","R2",NULL);
-
-  mb = cb->addMethod("union", &GModelFactory::booleanUnion);
-  mb->setDescription("returns a GModel that is the boolean union of 2 GModels. Third parameter tells if a new model has to be created");
-  mb->setArgNames("g1","g2","createNewGModel",NULL);
-
-  mb = cb->addMethod("intersection", &GModelFactory::booleanIntersection);
-  mb->setDescription("returns a GModel that is the boolean intersection of 2 GModels. Third parameter tells if a new model has to be created");
-  mb->setArgNames("g1","g2","createNewGModel",NULL);
-
-  mb = cb->addMethod("difference", &GModelFactory::booleanDifference);
-  mb->setDescription("returns a GModel that is the boolean difference of 2 GModels. Third parameter tells if a new model has to be created");
-  mb->setArgNames("g1","g2","createNewGModel",NULL);
-#if defined(HAVE_OCC)
-  cb = b->addClass<OCCFactory>("OCCFactory");
-  cb->setDescription("a GEntity Factory for OpenCascade");
-  mb = cb->setConstructor<OCCFactory>();
-  mb->setDescription("a GEntity Factory for OpenCascade");
-  cb->setParentClass<GModelFactory>();
-#endif
-}
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index 8e3b7ce6a3..218ae2cc87 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -1,7 +1,12 @@
+// 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 _GMODEL_FACTORY_H_
 #define _GMODEL_FACTORY_H_
-#include <vector>
 
+#include <vector>
 #include "GmshConfig.h"
 #include "fullMatrix.h"
 #include "SPoint3.h"
@@ -12,120 +17,90 @@ class GEdge;
 class GModel;
 class binding;
 
+// Abstract CAD creation factory.
 class GModelFactory {
- protected :
  public:
-  // constructor takes the GModel as input
   GModelFactory (){}
-  // creation methods
+  virtual ~GModelFactory(){}
+  // brep primitives
   enum arcCreationMethod {THREE_POINTS=1, CENTER_START_END=2};
   enum splineType {BEZIER=1, CATMULL_ROM=2};
-
-  virtual GVertex* createVertex (GModel *gm, double x, double y, double z, double lc) = 0;
-
-  virtual GEdge *createLine (GModel *, GVertex *v1, GVertex *v2) = 0;
-
-
-  virtual GEdge *createCircleArc (GModel *gm, const arcCreationMethod &method,
-				  GVertex *start, 
-				  GVertex *end, 
-				  const SPoint3 &aPoint)= 0;
-  inline GEdge *createCircleArcCenter (GModel *gm, double x, double y, double z,
-				       GVertex *start, 
-				       GVertex *end){
-    return createCircleArc (gm, CENTER_START_END,start,end,SPoint3(x,y,z));
-  }
-
-  inline GEdge *createCircleArc3Points (GModel *gm, double x, double y, double z,
-					GVertex *start, 
-					GVertex *end){
-    return createCircleArc (gm, THREE_POINTS,start,end,SPoint3(x,y,z));
-  }
-  virtual GEdge *createSpline (GModel *gm,const splineType &type,
-			       GVertex *start, 
-			       GVertex *end, 
-			       fullMatrix<double> *controlPoints)= 0;
-  inline GEdge *createBezier ( GModel *gm, GVertex *start, 
-			       GVertex *end, 
-			       fullMatrix<double> *controlPoints){
-    return createSpline(gm, BEZIER,start, end, controlPoints);
-  }
-
-  virtual GEntity* revolve (GModel *gm, GEntity*, 
-			    double x1, double y1, double z1, 
-			    double x2, double y2, double z2,
-			    double angle ) = 0;
-  virtual GEntity* extrude (GModel *gm, GEntity*, 
-			    double x1, double y1, double z1, 
-			    double x2, double y2, double z2) = 0;
-
-  inline GEntity* revolve_ (GModel *gm, GEntity* e, 
-			    double angle,
-			    fullMatrix<double> *axis){
-    revolve (gm, e, (*axis)(0,0),(*axis)(0,1),(*axis)(0,2),
-	     (*axis)(1,0),(*axis)(1,1),(*axis)(1,2),angle);
-  }
-  inline GEntity* extrude_ (GModel *gm, GEntity* e, 
-			    fullMatrix<double> *axis){
-    extrude (gm, e, (*axis)(0,0),(*axis)(0,1),(*axis)(0,2),
-	     (*axis)(1,0),(*axis)(1,1),(*axis)(1,2));
-  }
-  // primitives
-  virtual GEntity * sphere   (GModel *gm, double cx, double cy, double cz, double radius) = 0; 
-  virtual GEntity * cylinder (GModel *gm, std::vector<double> p1, std::vector<double> p2, 
-			      double radius) = 0; 
-  virtual GEntity * torus   (GModel *gm, std::vector<double> p1, std::vector<double> p2, 
-			     double radius1, double radius2) = 0; 
-  virtual GEntity * block  (GModel *gm, std::vector<double> p1, std::vector<double> p2) = 0;
-  virtual GEntity * cone   (GModel *gm, std::vector<double> p1, std::vector<double> p2,
-			    double radius1, double radius2) = 0; 
-
-  // boolean operators acting on 2 GEntities
-  virtual GModel * booleanUnion (GModel *obj, GModel*tool, int createNewModel) = 0;
-  virtual GModel * booleanIntersection (GModel *obj, GModel*tool, int createNewModel) = 0;
-  virtual GModel * booleanDifference (GModel *obj, GModel*tool, int createNewModel) = 0;
-  
-
-  // bindings
-  static void registerBindings(binding *b);
+  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 GEntity *revolve(GModel *gm, GEntity*, double x1, double y1, double z1, 
+                           double x2, double y2, double z2, double angle) = 0;
+  virtual GEntity *extrude(GModel *gm, GEntity*, double x1, double y1, double z1, 
+                           double x2, double y2, double z2) = 0;
+
+  // solid primitives
+  virtual GEntity *addSphere(GModel *gm, double cx, double cy, double cz, 
+                             double radius) = 0; 
+  virtual GEntity *addCylinder(GModel *gm, std::vector<double> p1, 
+                               std::vector<double> p2, double radius) = 0; 
+  virtual GEntity *addTorus(GModel *gm, std::vector<double> p1, 
+                            std::vector<double> p2, double radius1, 
+                            double radius2) = 0; 
+  virtual GEntity *addBlock(GModel *gm, std::vector<double> p1,
+                            std::vector<double> p2) = 0;
+  virtual GEntity *addCone(GModel *gm, std::vector<double> p1,
+                           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;
 };
 
 #if defined(HAVE_OCC)
+
 class OCCFactory : public GModelFactory {
  public:
-  OCCFactory ();
-  GVertex *createVertex (GModel *gm,double x, double y, double z, double lc);
-  virtual GEdge *createLine (GModel *gm,GVertex *v1, GVertex *v2);
-  GEdge *createCircleArc (GModel *gm,const arcCreationMethod &method,				  
-			  GVertex *start, 
-			  GVertex *end, 
-			  const SPoint3 &aPoint);
-  GEdge *createSpline (GModel *gm,const splineType &type,
-		       GVertex *start, 
-		       GVertex *end, 
-		       fullMatrix<double> *controlPoints);
-
-  GEntity* revolve (GModel *gm,GEntity*, 
-		    double x1, double y1, double z1, 
-		    double x2, double y2, double z2,
-		    double angle );
-
-  GEntity* extrude (GModel *gm,GEntity*, 
-		    double x1, double y1, double z1, 
-		    double x2, double y2, double z2);
-
-  GEntity * sphere   (GModel *gm,double cx, double cy, double cz, double radius); 
-  GEntity * cylinder (GModel *gm,std::vector<double> p1, std::vector<double> p2, 
-		      double radius); 
-  GEntity * torus   (GModel *gm,std::vector<double> p1, std::vector<double> p2, 
-		     double radius1, double radius2); 
-  GEntity * block  (GModel *gm,std::vector<double> p1, std::vector<double> p2); 
-  GEntity * cone   (GModel *gm,std::vector<double> p1, std::vector<double> p2, double radius1, double radius2); 
-  // booleans
-  GModel * booleanUnion (GModel *obj, GModel*tool, int createNewModel);
-  GModel * booleanIntersection (GModel *obj, GModel*tool, int createNewModel);
-  GModel * booleanDifference (GModel *obj, GModel*tool, int createNewModel);
+  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*, double x1, double y1, double z1, 
+                           double x2, double y2, double z2, double angle);
+  virtual GEntity *extrude(GModel *gm, GEntity*, double x1, double y1, double z1, 
+                           double x2, double y2, double z2);
+  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);
 };
+
 #endif
 
 #endif
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 59017b93cb..4ff8a85ad5 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -282,18 +282,26 @@ void OCC_Internals::healGeometry(double tolerance, bool fixsmalledges,
     
     if(sfwf->FixSmallEdges()){
       Msg::Info("- fixing wire frames");
-      if(sfwf->StatusSmallEdges(ShapeExtend_OK)) Msg::Info("no small edges found");
-      if(sfwf->StatusSmallEdges(ShapeExtend_DONE1)) Msg::Info("some small edges fixed");
-      if(sfwf->StatusSmallEdges(ShapeExtend_FAIL1)) Msg::Info("failed to fix some small edges");
+      if(sfwf->StatusSmallEdges(ShapeExtend_OK)) 
+        Msg::Info("no small edges found");
+      if(sfwf->StatusSmallEdges(ShapeExtend_DONE1))
+        Msg::Info("some small edges fixed");
+      if(sfwf->StatusSmallEdges(ShapeExtend_FAIL1)) 
+        Msg::Info("failed to fix some small edges");
     }
   
     if(sfwf->FixWireGaps()){
       Msg::Info("- fixing wire gaps");
-      if(sfwf->StatusWireGaps(ShapeExtend_OK)) Msg::Info("no gaps found");
-      if(sfwf->StatusWireGaps(ShapeExtend_DONE1)) Msg::Info("some 2D gaps fixed");
-      if(sfwf->StatusWireGaps(ShapeExtend_DONE2)) Msg::Info("some 3D gaps fixed");
-      if(sfwf->StatusWireGaps(ShapeExtend_FAIL1)) Msg::Info("failed to fix some 2D gaps");
-      if(sfwf->StatusWireGaps(ShapeExtend_FAIL2)) Msg::Info("failed to fix some 3D gaps");
+      if(sfwf->StatusWireGaps(ShapeExtend_OK))
+        Msg::Info("no gaps found");
+      if(sfwf->StatusWireGaps(ShapeExtend_DONE1))
+        Msg::Info("some 2D gaps fixed");
+      if(sfwf->StatusWireGaps(ShapeExtend_DONE2))
+        Msg::Info("some 3D gaps fixed");
+      if(sfwf->StatusWireGaps(ShapeExtend_FAIL1))
+        Msg::Info("failed to fix some 2D gaps");
+      if(sfwf->StatusWireGaps(ShapeExtend_FAIL2))
+        Msg::Info("failed to fix some 3D gaps");
     }
     
     shape = sfwf->Shape();
@@ -421,9 +429,9 @@ void OCC_Internals::loadSTEP(const char *fn)
 void OCC_Internals::writeSTEP(const char *fn)
 {
   STEPControl_Writer writer;
-  IFSelect_ReturnStatus status = writer.Transfer( shape, STEPControl_ManifoldSolidBrep ) ;
-  if ( status == IFSelect_RetDone ) 
-    status = writer.Write( (char*) fn) ;
+  IFSelect_ReturnStatus status = writer.Transfer(shape, STEPControl_ManifoldSolidBrep);
+  if (status == IFSelect_RetDone) 
+    status = writer.Write((char*)fn);
 }
 
 void OCC_Internals::loadIGES(const char *fn)
@@ -476,7 +484,7 @@ GEdge *OCC_Internals::addEdgeToModel(GModel *model, TopoDS_Edge edge)
 }
 
 GEdge *OCC_Internals::addEdgeToModel(GModel *model, TopoDS_Edge edge, GVertex *g1, 
-                                      GVertex *g2)
+                                     GVertex *g2)
 {
   OCCEdge *e = new OCCEdge(model, edge, model->maxEdgeNum() + 1, g1, g2);
   e->replaceEndingPoints (g1,g2);
@@ -815,7 +823,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 9bc4bc89bf..168360d4d0 100644
--- a/Geo/GModelIO_OCC.h
+++ b/Geo/GModelIO_OCC.h
@@ -48,7 +48,7 @@ class OCC_Internals {
   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);
+  void fillet(std::vector<TopoDS_Edge> &shapes, double radius);
   void applyBooleanOperator(TopoDS_Shape tool, const BooleanOperator &op);
 };
 
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 46d04e8e26..0ae4212d39 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -97,6 +97,8 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
 #include "Bindings.h"
 void linearSystemPETScRegisterBindings(binding *b) 
 {
+ // FIXME on complex arithmetic this crashes
+#if !defined(PETSC_USE_COMPLEX)
   classBinding *cb;
   methodBinding *cm;
   cb = b->addClass<linearSystemPETSc<PetscScalar> >("linearSystemPETSc");
@@ -107,7 +109,6 @@ void linearSystemPETScRegisterBindings(binding *b)
   cm->setArgNames(NULL);
   cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<PetscScalar> >::systemSolve);
   cm->setDescription("compute x = A^{-1}b");
-
   cb = b->addClass<linearSystemPETSc<fullMatrix<PetscScalar> > >("linearSystemPETScBlock");
   cb->setDescription("A linear system solver, based on PETSc");
   cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> >, int>();
@@ -115,6 +116,8 @@ void linearSystemPETScRegisterBindings(binding *b)
   cm->setArgNames("blockSize", NULL);
   cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<PetscScalar> >::systemSolve);
   cm->setDescription("compute x = A^{-1}b");
+#endif // FIXME
+
 }
 
 #endif // HAVE_PETSC
diff --git a/benchmarks/boolean/square1.lua b/benchmarks/boolean/square1.lua
new file mode 100644
index 0000000000..3fce43deaf
--- /dev/null
+++ b/benchmarks/boolean/square1.lua
@@ -0,0 +1,11 @@
+
+g = GModel()
+v1 = g:addVertex(0, 0, 0, 1)
+v2 = g:addVertex(1, 0, 0, 1)
+v3 = g:addVertex(1, 1, 0, 1)
+v4 = g:addVertex(0, 1, 0, 1)
+e1 = g:addLine(v1, v2)
+e2 = g:addLine(v2, v3)
+e3 = g:addLine(v3, v4)
+e4 = g:addLine(v4, v1)
+
diff --git a/benchmarks/boolean/square2.lua b/benchmarks/boolean/square2.lua
new file mode 100644
index 0000000000..ff0236efcd
--- /dev/null
+++ b/benchmarks/boolean/square2.lua
@@ -0,0 +1,15 @@
+
+g = GModel()
+v1 = g:addVertex(0, 0, 0, 1)
+v2 = g:addVertex(1, 0, 0, 1)
+e1 = g:addLine(v1, v2)
+
+dir = fullMatrix(2,3)
+dir:set(0,0, 0);
+dir:set(0,1, 0);
+dir:set(0,2, 0);
+dir:set(1,0, 0);
+dir:set(1,1, 1);
+dir:set(1,2, 0);
+f1 = g:extrude(e1, dir)
+
diff --git a/benchmarks/boolean/wikipedia.lua b/benchmarks/boolean/wikipedia.lua
new file mode 100644
index 0000000000..b6494b4320
--- /dev/null
+++ b/benchmarks/boolean/wikipedia.lua
@@ -0,0 +1,26 @@
+
+-- from http://en.wikipedia.org/wiki/Constructive_solid_geometry
+
+R = 1.0;
+s = .8;
+t = 1.35;
+myModel = GModel();
+myModel:addBlock({-R,-R,-R},{R,R,R});
+
+myTool = GModel();
+myTool:addSphere(0,0,0,R*t);
+
+myModel:computeIntersection(myTool,0);
+
+myTool2 = GModel();
+myTool2:addCylinder({-2*R,0,0},{2*R,0,0},R*s);
+
+myTool3 = GModel();
+myTool3:addCylinder({0,-2*R,0},{0,2*R,0},R*s);
+
+myModel2 = GModel();
+myModel2:addCylinder({0,0,-2*R},{0,0,2*R},R*s);
+myModel2:computeUnion(myTool2,0);
+myModel2:computeUnion(myTool3,0);
+
+myModel:computeDifference(myModel2,0);
-- 
GitLab