diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index d8c773c4f8b7eae18272cf7b1c070907c1082187..85a0163ea8dd4883eca77701f8de32d63ae0ece7 100644
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -144,16 +144,12 @@ void Cell::findBdElement(int i, std::vector<MVertex*>& vertices) const
         vertices.push_back(_v[MTetrahedron::faces_tetra(i, j)]);
       return;
     case 5:
-      if(i < 3) {
-      for(int j = 0; j < 3; j++)
-        vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
-      }
-      else {
-        vertices.push_back(_v[0]);
-        vertices.push_back(_v[3]);
-        vertices.push_back(_v[2]);
-        vertices.push_back(_v[1]);
-      }
+      if(i < 4)
+        for(int j = 0; j < 3; j++)
+          vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
+      else
+        for(int j = 0; j < 4; j++)
+          vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
       return;
     case 6:
       if(i < 2)
@@ -270,8 +266,145 @@ int Cell::findBdCellOrientation(Cell* cell, int i) const
          _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[1]->getNum() &&
          _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[0]->getNum())
         return -1;
-    case 5: return 0;
-    case 6: return 0;
+    case 5:
+      if(i < 4) {
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum())
+          return 1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum())
+          return 1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum())
+          return 1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum())
+          return -1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum())
+          return -1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum())
+          return -1;
+      }
+      else {
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[3]->getNum())
+          return 1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[3]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[0]->getNum())
+          return 1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[3]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[1]->getNum())
+          return 1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[3]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[2]->getNum())
+          return 1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[3]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[1]->getNum())
+          return -1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[3]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[0]->getNum())
+          return -1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[3]->getNum())
+          return -1;
+        if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[3]->getNum() &&
+           _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[2]->getNum())
+          return -1;
+      }
+      return 0;
+    case 6:
+      if(i < 2) {
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum())
+          return 1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum())
+          return 1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum())
+          return 1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum())
+          return -1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum())
+          return -1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum())
+          return -1;
+      }
+      else {
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[3]->getNum())
+          return 1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[3]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[0]->getNum())
+          return 1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[3]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[1]->getNum())
+          return 1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[3]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[2]->getNum())
+          return 1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[3]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[1]->getNum())
+          return -1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[3]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[0]->getNum())
+          return -1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[3]->getNum())
+          return -1;
+        if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
+           _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
+           _v[MPrism::faces_prism(i, 2)]->getNum() == v[3]->getNum() &&
+           _v[MPrism::faces_prism(i, 3)]->getNum() == v[2]->getNum())
+          return -1;
+      }
     case 8:
       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() &&
          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() &&
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 6638f0fd7495df34a7ae33ec593adb73c4951a53..85990e92fc31c9d47eb1bd9c05e0a8735892ac5c 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -82,11 +82,12 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
     MElement* element = elements.at(i);
     int dim = element->getDim();
     int type = element->getType();
-    if(type == TYPE_PYR || type == TYPE_PRI ||
+    if(//type == TYPE_PYR || type == TYPE_PRI ||
        type == TYPE_POLYG || type == TYPE_POLYH) {
       Msg::Error("Mesh element type %d not implemented in homology solver", type);
     }
-    if(type == TYPE_QUA || type == TYPE_HEX)
+    if(type == TYPE_QUA || type == TYPE_HEX ||
+       type == TYPE_PYR || type == TYPE_PRI)
       _simplicial = false;
     std::pair<Cell*, bool> maybeCell = Cell::createCell(element, domain);
     if(!maybeCell.second) {
diff --git a/Geo/Chain.cpp b/Geo/Chain.cpp
index 3c7b04b3f6f74f972b6077de4758a7f2c26cef2d..13ba5623c33bb14b1402d43a5aa6c570083a2adc 100644
--- a/Geo/Chain.cpp
+++ b/Geo/Chain.cpp
@@ -222,16 +222,12 @@ void ElemChain::getBoundaryVertices(int i, int dim, int numVertices,
         vertices.push_back(v[MTetrahedron::faces_tetra(i, j)]);
       return;
     case 5:
-      if(i < 3) {
-      for(int j = 0; j < 3; j++)
-        vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
-      }
-      else {
-        vertices.push_back(v[0]);
-        vertices.push_back(v[3]);
-        vertices.push_back(v[2]);
-        vertices.push_back(v[1]);
-      }
+      if(i < 3)
+        for(int j = 0; j < 3; j++)
+          vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
+      else
+        for(int j = 0; j < 4; j++)
+          vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
       return;
     case 6:
       if(i < 2)
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 114f90766ee80676893d6a313a35048c333e2348..d0da0319e7d68d149620ba9e4a98d80584db5921 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -337,7 +337,7 @@ void ChainComplex::computeHomology(bool dual)
     }
 
     // 2) this dimension is empty
-    else if(getHMatrix(lowDim) == NULL){
+    else if(getHMatrix(setDim) == NULL){
       setHbasis(setDim, NULL);
     }
     // 3) No higher dimension cells -> none of the cycles are boundaries
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index 650cb5115ba95f254e02a574324223988d9ba080..bd553a6e60f9b557a84c41f6cbc5cd1d3be7edc8 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -27,6 +27,12 @@ int MPyramid::getVolumeSign()
   else return 0;
 }
 
+void MPyramid::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+  *npts = getNGQPyrPts(pOrder);
+  *pts = getGQPyrPts(pOrder);
+}
+
 const nodalBasis* MPyramid::getFunctionSpace(int o) const
 {
   int order = (o == -1) ? getPolynomialOrder() : o;
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index bad11f2bea2b9d22143221a31ddaffb217f05c85..2403d674ea74acdef3776761766ad27fac476657 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -149,6 +149,7 @@ class MPyramid : public MElement {
       return false;
     return true;
   }
+  void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
   static int edges_pyramid(const int edge, const int vert)
   {
     static const int e[8][2] = {