diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index f08f804863698311127a32b1e386b97ba26ebe7c..43fe44cdd134200fa898cb1d32156800187001dc 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -633,59 +633,11 @@ int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
   return count;
 }
 
-std::vector<int> CellComplex::bettiCoreduceComplex()
+void CellComplex::bettiReduceComplex()
 {
-  Msg::Debug("Cell Complex betti coreduction:");
-  Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
-            getSize(3), getSize(2), getSize(1), getSize(0));
-
-  std::vector<int> betti(4,0);
-  if(!getSize(0)) return betti;
-
-  std::vector<Cell*> empty;
-  if(relative()) {
-    removeSubdomain();
-    int count = 0;
-    for(int dim = 0; dim < 4; dim++){
-      citer cit = firstCell(dim);
-      while(cit != lastCell(dim)){
-        Cell* cell = *cit;
-        int count =+ coreduction(cell, -1, empty);
-        if(count != 0) break;
-        cit++;
-      }
-    }
-    for(int j = 1; j <= getDim(); j++) count += coreduction(j, -1, empty);
-  }
-
-  for(int i = 0; i < 4; i++) {
-    while (getSize(i) != 0){
-      citer cit = firstCell(i);
-      Cell* cell = *cit;
-
-      Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
-                 getSize(3), getSize(2), getSize(1), getSize(0));
-
-      int count = 1;
-
-      removeCell(cell, false);
-      betti.at(i)++;
-
-      count += coreduction(cell, -1, empty);
-      for(int j = 1; j <= getDim(); j++) count += coreduction(j, -1, empty);
-
-      std::string domainstr = "";
-      getDomain(cell, domainstr);
-
-      Msg::Debug("Omitted %d-cell in %s that caused %d reductions",
-                 cell->getDim(), domainstr.c_str(), count);
-      Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
-                 getSize(3), getSize(2), getSize(1), getSize(0));
-    }
-  }
-
-  _reduced = true;
-  return betti;
+  reduceComplex(3, true);
+  for(int i = 1; i <= 3; i++) cocombine(i-1);
+  return;
 }
 
 int CellComplex::combine(int dim)
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index e80238ce42844496fd00ad11815e27d18edc8056..830e7f4767e31cfa52ab0cd71a1d5ddb52ad3c63 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -170,8 +170,8 @@ class CellComplex
   int reduceComplex(int combine=1, bool omit=true, bool homseq=false);
   int coreduceComplex(int combine=1, bool omit=true, int heuristic=0);
 
-  // compute Betti numbers of the cell complex using coreduction
-  std::vector<int> bettiCoreduceComplex();
+  // reduce cell complex for Betti number computation
+  void bettiReduceComplex();
 
   bool isReduced() const { return _reduced; }
 
diff --git a/Geo/Chain.h b/Geo/Chain.h
index bc3d91073a0b9b3013b4618b6ac6d2f12cfa228a..1a0086afaff8030c9b0088ff9359dd49bd68baac 100644
--- a/Geo/Chain.h
+++ b/Geo/Chain.h
@@ -257,6 +257,7 @@ public:
 template <class C>
 Chain<C>::Chain(GModel* m, int physicalGroup)
 {
+  _dim = 0;
   std::vector<int> groups(1, physicalGroup);
   std::vector<GEntity*> entities;
   findEntitiesInPhysicalGroups(m, groups, entities);
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 3a52387cf6b267a74d75b197a582460cd3aa036f..2733322dd905df60f3068a9de991e02d759a9bd5 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -519,17 +519,27 @@ void Homology::findBettiNumbers()
     if(_cellComplex == NULL) _createCellComplex();
     if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
 
-    Msg::StatusBar(true, "Computing betti numbers...");
-
+    Msg::StatusBar(true, "Reducing cell complex...");
     double t1 = Cpu();
 
-    std::vector<int> betti = _cellComplex->bettiCoreduceComplex();
+    _cellComplex->bettiReduceComplex();
 
     double t2 = Cpu();
+    Msg::StatusBar(true, "Done reducing cell complex (%g s)", t2 - t1);
+    Msg::Info("%d volumes, %d faces, %d edges, and %d vertices",
+              _cellComplex->getSize(3), _cellComplex->getSize(2),
+              _cellComplex->getSize(1), _cellComplex->getSize(0));
 
-    Msg::StatusBar(true, "Betti numbers computed (%g s)", t2 - t1);
 
-    for(int i = 0; i < 4; i++) _betti[i] = betti.at(i);
+    Msg::StatusBar(true, "Computing betti numbers...");
+    t1 = Cpu();
+    ChainComplex chainComplex = ChainComplex(_cellComplex);
+    chainComplex.computeHomology();
+
+    for(int i = 0; i < 4; i++) _betti[i] = chainComplex.getBasisSize(i, 3);
+
+    t2 = Cpu();
+    Msg::StatusBar(true, "Betti numbers computed (%g s)", t2 - t1);
   }
 
   std::string domain = _getDomainString(_domain, _subdomain);