diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 0e19ec7177682627bc5894399ba75c0f3643d79f..e4345d5c65cf51e0fb238ec49c03dcb7e9ad205a 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -2214,12 +2214,9 @@ bool OCC_Internals::booleanOperator(int tag, BooleanOperator op,
     return false;
   }
 
-  // don't try to preserve numbering if we specify the tag explicitly, or if
-  // there is a problem
-  bool bug1 = (objectDimTags.size() + toolDimTags.size() != mapModified.size());
-  bool bug2 = (op == OCC_Internals::Union); // strange fuse behavior in OCC 7.1
-  if(tag >= 0 || bug1 || bug2){
-    if(bug1) Msg::Error("Wrong shape count in boolean operation");
+  // if we specify the tag explicitly, just go ahead and bind the resulting
+  // shape (and sub-shapes)
+  if(tag >= 0){
     if(removeObject){
       for(unsigned int i = 0; i < objectDimTags.size(); i++){
         int d = objectDimTags[i].first;
@@ -2239,26 +2236,22 @@ bool OCC_Internals::booleanOperator(int tag, BooleanOperator op,
     return true;
   }
 
-  // otherwise, try to preserve the numbering
+  // otherwise, preserve the numbering of the input shapes that did not change,
+  // or that were replaced by a single shape
   for(unsigned int i = 0; i < objectDimTags.size(); i++){
     int dim = objectDimTags[i].first;
     int tag = objectDimTags[i].second;
-    if(mapDeleted[i] && !mapGenerated[i].Extent()){
-      // the shape has been deleted
-      if(removeObject && _isBound(dim, mapOriginal[i])){
-        unbind(mapOriginal[i], dim, tag, true);
-      }
+    if(mapDeleted[i]){ // object deleted
+      if(removeObject) unbind(mapOriginal[i], dim, tag, true);
     }
-    else if(mapModified[i].Extent() == 0){
-      // the shape has not been modified
+    else if(mapModified[i].Extent() == 0){ // object not modified
       outDimTags.push_back(std::pair<int, int>(dim, tag));
+      unbind(mapOriginal[i], dim, tag, false); // not recursive!
+      bind(mapOriginal[i], dim, tag, false); // not recursive!
     }
-    else if(mapModified[i].Extent() == 1){
+    else if(mapModified[i].Extent() == 1){ // object replaced by single one
       if(removeObject){
-        // the shape has been replaced by a single shape, keep the same tag
-        if(_isBound(dim, mapOriginal[i])){
-          unbind(mapOriginal[i], dim, tag, true);
-        }
+        unbind(mapOriginal[i], dim, tag, true);
         bind(mapModified[i].First(), dim, tag, false); // not recursive!
         int t = _find(dim, mapModified[i].First());
         if(tag != t)
@@ -2267,32 +2260,24 @@ bool OCC_Internals::booleanOperator(int tag, BooleanOperator op,
       }
     }
     else{
-      if(removeObject && _isBound(dim, mapOriginal[i])){
-        unbind(mapOriginal[i], dim, tag, true);
-      }
+      if(removeObject) unbind(mapOriginal[i], dim, tag, true);
     }
   }
-
   for(unsigned int i = 0; i < toolDimTags.size(); i++){
     int k = objectDimTags.size() + i;
     int dim = toolDimTags[i].first;
     int tag = toolDimTags[i].second;
-    if(mapDeleted[k] && !mapGenerated[k].Extent()){
-      // the shape has been deleted
-      if(removeTool && _isBound(dim, mapOriginal[k])){
-        unbind(mapOriginal[k], dim, tag, true);
-      }
+    if(mapDeleted[k]){ // tool deleted
+      if(removeTool) unbind(mapOriginal[k], dim, tag, true);
     }
-    else if(mapModified[k].Extent() == 0){
-      // the shape has not been modified
+    else if(mapModified[k].Extent() == 0){ // tool not modified
       outDimTags.push_back(std::pair<int, int>(dim, tag));
+      unbind(mapOriginal[k], dim, tag, false); // not recursive!
+      bind(mapOriginal[k], dim, tag, false); // not recursive!
     }
     else if(mapModified[k].Extent() == 1){
-      if(removeTool){
-        // the shape has been replaced by a single shape, keep the same tag
-        if(_isBound(dim, mapOriginal[k])){
-          unbind(mapOriginal[k], dim, tag, true);
-        }
+      if(removeTool){ // tool replaced by single one
+        unbind(mapOriginal[k], dim, tag, true);
         bind(mapModified[k].First(), dim, tag, false); // not recursive!
         int t = _find(dim, mapModified[k].First());
         if(tag != t)
@@ -2301,14 +2286,12 @@ bool OCC_Internals::booleanOperator(int tag, BooleanOperator op,
       }
     }
     else{
-      if(removeTool && _isBound(dim, mapOriginal[k])){
-        unbind(mapOriginal[k], dim, tag, true);
-      }
+      if(removeTool) unbind(mapOriginal[k], dim, tag, true);
     }
   }
 
-  // bind all remaining entities and add (only) new one to the returned list
-  _multiBind(result, tag, outDimTags, false, true, true);
+  // bind all remaining entities and add the new ones to the returned list
+  _multiBind(result, -1, outDimTags, false, true, true);
   _filterTags(outDimTags, minDim);
 
   return true;
@@ -2383,6 +2366,8 @@ bool OCC_Internals::_transform(const std::vector<std::pair<int, int> > &inDimTag
       }
       result = gtfo->Shape();
     }
+    // FIXME we should implement rebind(object, result, dim) which would
+    // unbind/bind all subshapes to the same tags
     unbind(object, dim, tag, true);
     bind(result, dim, tag, true);
   }
diff --git a/benchmarks/3d/Cube-05.geo b/benchmarks/3d/Cube-05.geo
index 12e4015f450909bc919860a9191fd9346cc00ce0..720725b92804e7648635f4a553a7bb5f6237e472 100644
--- a/benchmarks/3d/Cube-05.geo
+++ b/benchmarks/3d/Cube-05.geo
@@ -1,5 +1,5 @@
-lv = .1;
-lc = .1;
+lv = .15;
+lc = .12;
 Point(1) = {0.0,0.0,0.0,lv};
 Point(2) = {1,0.0,0.0,lv};
 Point(3) = {1,1,0.0,lv};
diff --git a/demos/boolean/baffles.geo b/demos/boolean/baffles.geo
index f3a2c1948cd62f0caf47cc31da8d31bfb13cee66..2b156105966fc3647d1f54d935fdc304962a22f1 100644
--- a/demos/boolean/baffles.geo
+++ b/demos/boolean/baffles.geo
@@ -1,8 +1,8 @@
 SetFactory("OpenCASCADE");
 Block(1) = {-0, 0, 0, 1, 1, 1};
+p() = PointsOf{ Volume{1}; };
+Characteristic Length{p()} = 0.2;
 
-
-// intersecting:
 Rectangle(7) = {0.4, 0.7, 0.2, 0.3, 0.3, 0};
 Rectangle(8) = {-0.2, 0.3, 0.5, 0.3, 0.3, 0};
 Rectangle(9) = {0.8, -0.2, 0.5, 0.3, 0.3, 0};
@@ -11,6 +11,5 @@ Rectangle(10) = {0.3, 0.3, 0.5, 0.3, 0.3, 0};
 a() = BooleanFragments{ Volume{1}; Delete; }{ Surface{7:10}; Delete; };
 Printf("a = ", a());
 
-//Characteristic Length{Point "*"} = 0.1;
-//p() = PointsOf{ Surface{a({1:#a()-1})}; };
-//Characteristic Length{p()} = 0.01;
+p() = PointsOf{ Surface{a({1:6})}; };
+Characteristic Length{p()} = 0.04;