diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index cf39f61202a9d4d05ccbd625bb3df18c05eb7a02..29b31d0bc4e4f8ae3271af1c5a105442d1a52ea6 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -498,16 +498,17 @@ int CellComplex::reduceComplex(int combine, bool omit, bool homseq)
   if(combine > 0) this->combine(3);
 
   if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
-  else reduction(2, -1, empty);
+  else if(combine > 1) reduction(2, -1, empty);
 
   if(combine > 0) this->combine(2);
 
   if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
-  else reduction(1, -1, empty);
+  else if(combine > 1) reduction(1, -1, empty);
 
   if(combine > 0) this->combine(1);
 
   if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
+  else if(combine > 1) reduction(0, -1, empty);
 
   Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
              getSize(3), getSize(2), getSize(1), getSize(0));
@@ -604,17 +605,17 @@ int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
   if(combine > 0) this->cocombine(0);
 
   if(combine > 2) for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
-  else coreduction(1, -1, empty);
+  else if(combine > 1) coreduction(1, -1, empty);
 
   if(combine > 0) this->cocombine(1);
 
   if(combine > 2)  for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
-  else coreduction(2, -1, empty);
+  else if(combine > 1) coreduction(2, -1, empty);
 
   if(combine > 0) this->cocombine(2);
 
   if(combine > 2) for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
-  else coreduction(3, -1, empty);
+  else if(combine > 1) coreduction(3, -1, empty);
 
   coherent();
   Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
@@ -949,6 +950,16 @@ Cell* CellComplex::getACell(int dim, int domain)
 bool CellComplex::restoreComplex()
 {
   if(_saveorig){
+
+    for(unsigned int i = 0; i < _removedcells.size(); i++) {
+      Cell* cell = _removedcells.at(i);
+      if(cell->isCombined()) {
+        delete cell;
+        _deleteCount++;
+      }
+    }
+    _removedcells.clear();
+
     for(int i = 0; i < 4; i++){
 
       for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index eba01b7e33c7e2f38ab184a68b0aa21b3db615af..7b8375bbaad3d24e41bc33f3353d187f0641c14d 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -256,9 +256,10 @@ void ChainComplex::Inclusion(int lowDim, int highDim)
   destroy_gmp_normal_form(normalForm);
 }
 
-void ChainComplex::Quotient(int dim)
+void ChainComplex::Quotient(int dim, int setDim)
 {
   if(dim < 0 || dim > 4 || _JMatrix[dim] == NULL) return;
+  if(setDim < 0 || setDim > 4) return;
 
   gmp_matrix* JMatrix =
     copy_gmp_matrix(_JMatrix[dim], 1, 1,
@@ -283,10 +284,12 @@ void ChainComplex::Quotient(int dim)
       destroy_gmp_normal_form(normalForm);
       return;
     }
-    if(mpz_cmp_si(elem,1) > 0) _torsion[dim].push_back(mpz_get_si(elem));
+    if(mpz_cmp_si(elem,1) > 0) {
+      _torsion[setDim].push_back(mpz_get_si(elem));
+    }
   }
 
-  int rank = cols - _torsion[dim].size();
+  int rank = cols - _torsion[setDim].size();
   if(rows - rank > 0){
     gmp_matrix* Hbasis =
       copy_gmp_matrix(normalForm->left, 1, rank+1, rows, rows);
@@ -304,6 +307,8 @@ void ChainComplex::computeHomology(bool dual)
   int highDim = 0;
   int setDim = 0;
 
+  if(dual) transposeHMatrices();
+
   for(int i=-1; i < 4; i++){
 
     if(dual){
@@ -361,7 +366,7 @@ void ChainComplex::computeHomology(bool dual)
 		      create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) );
       }
       Inclusion(lowDim, highDim);
-      Quotient(lowDim);
+      Quotient(lowDim, setDim);
 
       if(getCodHMatrix(highDim) == NULL){
         setHbasis(setDim,
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 04ad96fb82323dc0ec3760bf37d9a89e969c7def..62fd8726ae33a045d9b13ddcb749fb48576ea9dc 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -129,7 +129,7 @@ class ChainComplex
   // compute quotient problem for given inclusion relation j to find
   // representatives of homology group generators and possible
   // torsion coeffcients
-  void Quotient(int dim);
+  void Quotient(int dim, int setDim);
 
   // transpose the boundary operator matrices, these are boundary operator
   // matrices for the dual mesh
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 1e8f66e340ba018cbdaa8bbae58ab18333642462..41c7716a6384ef6fb1e6aea15186fa5151a401db 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -290,7 +290,6 @@ void Homology::findCohomologyBasis(std::vector<int> dim)
   Msg::StatusBar(true, "Computing cohomology space bases ...");
   t1 = Cpu();
   ChainComplex chainComplex = ChainComplex(_cellComplex);
-  chainComplex.transposeHMatrices();
   chainComplex.computeHomology(true);
   t2 = Cpu();
   Msg::StatusBar(true, "Done computing cohomology space bases (%g s)", t2- t1);
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 78d0ba427c2e273c5b50e4583bfcc84654d6ae4b..0dd4a1d4b4afd23e374c4f97af4c8c5497947dfd 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -138,6 +138,9 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
       homology->addChainsToModel(dim, pviews, hompg);
       if(hompg != -1) hompg += homology->betti(dim);
     }
+  }
+  for(unsigned int i = 0; i < dimsave.size(); i++) {
+    int dim = dimsave.at(i);
     if(dim > -1 && dim < 4 && coh != 0) {
       homology->addCochainsToModel(dim, pviews, cohpg);
       if(cohpg != -1) cohpg += homology->betti(dim);