diff --git a/Geo/Chain.h b/Geo/Chain.h
index f8cc1ca9ca3bbf91fefee308651f38984729cd01..45fbee74dd127ed7471a8e8b4f232ce81c3fb600 100644
--- a/Geo/Chain.h
+++ b/Geo/Chain.h
@@ -171,35 +171,101 @@ template <class C>
 class Chain : public VectorSpaceCat<Chain<C>, C>
 {
 private:
+  // Dimension of the chain
   int _dim;
+  // Elementary chains and their coefficients in the chain
   std::map<ElemChain, C> _elemChains;
+  // A name for the chain
   std::string _name;
 
 public:
+  // Elementary chain iterators
   typedef typename std::map<ElemChain, C>::iterator ecit;
   typedef typename std::map<ElemChain, C>::const_iterator cecit;
 
+  // Create zero chain
   Chain() : _dim(-1), _name("") {}
+
+  // Create chain from Gmsh model physical group
+  // (all mesh elements in the physical group are treated as
+  //  elementary chains with coefficient 1)
+  Chain(GModel* m, int physicalGroup);
+
+  // Get/set the chain name
   std::string getName() const { return _name; }
   void setName(std::string name) { _name = name; }
+
+  // Get chain dimension
   int getDim() const { return _dim; }
+
+  // True if a zero element of a chain space
   bool isZero() const { return _elemChains.empty(); }
+
+  // Iterators to elementrary chains in the chain
   cecit firstElemChain() const { return _elemChains.begin(); }
   cecit lastElemChain() const { return _elemChains.end(); }
 
+  // Add mesh element or elementary chain with given coefficient to the chain
   void addMeshElement(MElement* e, C coeff=1);
   void addElemChain(const ElemChain& c, C coeff=1);
 
+  // Vector space operations for chains (these two induce the rest)
   Chain<C>& operator+=(const Chain<C>& chain);
   Chain<C>& operator*=(const C& coeff);
 
+  // Get elementary chain coefficient the chain
   C getCoefficient(const ElemChain& c2) const;
+
+  // Get mesh element (or its indicated face, edge, or vertex)
+  // coefficient in the chain, interpreted as a elementary chain
+  C getCoefficient(MElement* e, int subElement=-1) const;
+
+  // Get the boundary chain of this chain
   Chain<C> getBoundary() const;
-  Chain<C> getTrace(const std::vector<GEntity*>& subdomain) const;
 
+  // Get a chain which contains elementary chains that are
+  // fully in the given physical group
+  Chain<C> getTrace(GModel* m, int physicalGroup) const;
+  Chain<C> getTrace(GModel* m, const std::vector<int>& physicalGroups) const;
+
+  // Add chain to Gmsh model as a physical group,
+  // elementary chains are turned into mesh elements with
+  // orientation and multiplicity given by elementary chain coefficient
+  // (and create a post-processing view)
   void addToModel(GModel* m, bool post=true) const;
 };
 
+template <class C>
+Chain<C>::Chain(GModel* m, int physicalGroup)
+{
+  std::vector<GEntity*> entities;
+  std::map<int, std::vector<GEntity*> > groups[4];
+  m->getPhysicalGroups(groups);
+  std::map<int, std::vector<GEntity*> >::iterator it;
+
+  for(int j = 0; j < 4; j++){
+    it = groups[j].find(physicalGroup);
+    if(it != groups[j].end()){
+      _dim = j;
+      std::vector<GEntity*> physicalGroup = it->second;
+      for(unsigned int k = 0; k < physicalGroup.size(); k++){
+        entities.push_back(physicalGroup.at(k));
+      }
+    }
+  }
+  if(entities.empty()) {
+    Msg::Error("Physical group %d does not exist", physicalGroup);
+  }
+  for(unsigned int i = 0; i < entities.size(); i++) {
+    GEntity* e = entities.at(i);
+    for(unsigned int j = 0; j < e->getNumMeshElements(); j++) {
+      this->addMeshElement(e->getMeshElement(j));
+    }
+    this->setName(m->getPhysicalName(this->getDim(),
+                                     physicalGroup));
+  }
+}
+
 template <class C>
 C Chain<C>::getCoefficient(const ElemChain& c2) const
 {
@@ -209,6 +275,37 @@ C Chain<C>::getCoefficient(const ElemChain& c2) const
   else return 0;
 }
 
+template <class C>
+C Chain<C>::getCoefficient(MElement* e, int subElement) const
+{
+  if(this->getDim() == e->getDim()) {
+    ElemChain ec(e);
+    return this->getCoefficient(ec);
+  }
+  if(subElement == -1) return 0;
+  std::vector<MVertex*> v;
+  if(this->getDim() == 0) {
+    if(subElement >= e->getNumVertices()) return 0;
+    v = std::vector<MVertex*>(1, e->getVertex(subElement));
+  }
+  else if(this->getDim() == 1) {
+    if(subElement >= e->getNumEdges()) return 0;
+    e->getEdgeVertices(subElement, v);
+    v.resize(2);
+  }
+  else if(this->getDim() == 2) {
+    if(subElement >= e->getNumFaces()) return 0;
+    e->getFaceVertices(subElement, v);
+    if(e->getType() == TYPE_TET ||
+       (e->getType() == TYPE_PRI && subElement < 4) ||
+       (e->getType() == TYPE_PYR && subElement < 2))
+      v.resize(3);
+    else v.resize(4);
+  }
+  ElemChain ec(this->getDim(), v);
+  return this->getCoefficient(ec);
+}
+
 template <class C>
 C incidence(const Chain<C>& c1, const Chain<C>& c2)
 {
@@ -252,20 +349,50 @@ Chain<C> Chain<C>::getBoundary() const
 }
 
 template <class C>
-Chain<C> Chain<C>::getTrace(const std::vector<GEntity*>& subdomain) const
+Chain<C> Chain<C>::getTrace(GModel* m, int physicalGroup) const
+{
+  std::vector<int> groups(1, physicalGroup);
+  return this->getTrace(m, groups);
+}
+
+template <class C>
+Chain<C> Chain<C>::getTrace(GModel* m,
+                            const std::vector<int>& physicalGroups) const
 {
   Chain<C> result;
-  for(cecit it = _elemChains.begin(); it != _elemChains.end(); it++) {
+  std::map<int, std::vector<GEntity*> > groups[4];
+  m->getPhysicalGroups(groups);
+  std::map<int, std::vector<GEntity*> >::iterator it;
+  std::vector<GEntity*> entities;
+  for(unsigned int i = 0; i < physicalGroups.size(); i++){
+    bool found = false;
+    for(int j = 0; j < 4; j++){
+      it = groups[j].find(physicalGroups.at(i));
+      if(it != groups[j].end()){
+        found = true;
+        std::vector<GEntity*> physicalGroup = it->second;
+        for(unsigned int k = 0; k < physicalGroup.size(); k++){
+          entities.push_back(physicalGroup.at(k));
+        }
+      }
+    }
+    if(!found) {
+      Msg::Error("Physical group %d does not exist",
+                 physicalGroups.at(i));
+    }
+    if(entities.empty()) return result;
+    for(cecit it = _elemChains.begin(); it != _elemChains.end(); it++) {
     bool inDomain = false;
-    for(unsigned int i = 0; i < subdomain.size(); i++) {
-      if(it->first.inEntity(subdomain.at(i))) {
+    for(unsigned int i = 0; i < entities.size(); i++) {
+      if(it->first.inEntity(entities.at(i))) {
         inDomain = true;
         break;
       }
     }
     if(inDomain) result.addElemChain(it->first, it->second);
+    }
+    return result;
   }
-  return result;
 }
 
 template <class C>
diff --git a/Plugin/HomologyPostProcessing.cpp b/Plugin/HomologyPostProcessing.cpp
index 3511f97881a41c948ce30ae744cd0eb082f28fc4..10bcf9ef91a9480fb8514568cd53e73d58a7cdd1 100644
--- a/Plugin/HomologyPostProcessing.cpp
+++ b/Plugin/HomologyPostProcessing.cpp
@@ -18,7 +18,6 @@
 #if defined(HAVE_KBIPACK)
 
 StringXNumber HomologyPostProcessingOptions_Number[] = {
-  {GMSH_FULLRC, "DimensionOfChains", NULL, 1},
   {GMSH_FULLRC, "ApplyBoundaryOperatorToResults", NULL, 0}
 };
 
@@ -92,23 +91,6 @@ bool GMSH_HomologyPostProcessingPlugin::parseStringOpt
   return true;
 }
 
-void GMSH_HomologyPostProcessingPlugin::createChains
-(int dim, GModel* m,
- const std::vector<GEntity*>& chainEntities,
- const std::vector<std::string>& chainNames,
- std::vector<Chain<int> >& chains)
-{
-  for(unsigned int i = 0; i < chainEntities.size(); i++) {
-    GEntity* e = chainEntities.at(i);
-    Chain<int> chain;
-    for(unsigned int j = 0; j < e->getNumMeshElements(); j++) {
-      chain.addMeshElement(e->getMeshElement(j));
-    }
-    chains.push_back(chain);
-    chains.at(i).setName(chainNames.at(i));
-  }
-}
-
 int GMSH_HomologyPostProcessingPlugin::detIntegerMatrix
 (std::vector<int>& matrix)
 {
@@ -145,13 +127,7 @@ PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
   std::string opString2 = HomologyPostProcessingOptions_String[2].def;
   std::string cname = HomologyPostProcessingOptions_String[4].def;
   std::string traceString = HomologyPostProcessingOptions_String[3].def;
-  int dim = (int)HomologyPostProcessingOptions_Number[0].def;
-  int bd = (int)HomologyPostProcessingOptions_Number[1].def;
-
-  if(dim < 0 || dim > 3) {
-    Msg::Error("Invalid cell dimension %d", dim);
-    return 0;
-  }
+  int bd = (int)HomologyPostProcessingOptions_Number[0].def;
 
   GModel* m = GModel::current();
 
@@ -227,52 +203,20 @@ PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
   std::vector<int> tracePhysicals;
   if(!parseStringOpt(3, tracePhysicals)) return 0;
 
-  std::vector<GEntity*> basisEntities;
-  std::vector<GEntity*> basisEntities2;
-  std::map<int, std::vector<GEntity*> > groups[4];
-  m->getPhysicalGroups(groups);
-  std::map<int, std::vector<GEntity*> >::iterator it;
-
-  std::vector<std::string> curBasisNames;
+  std::vector<Chain<int> > curBasis;
   for(unsigned int i = 0; i < basisPhysicals.size(); i++){
-    curBasisNames.push_back(m->getPhysicalName(dim, basisPhysicals.at(i)));
-    it = groups[dim].find(basisPhysicals.at(i));
-    if(it != groups[dim].end()){
-      std::vector<GEntity*> physicalGroup = it->second;
-      for(unsigned int k = 0; k < physicalGroup.size(); k++){
-        basisEntities.push_back(physicalGroup.at(k));
-        break;
-      }
-    }
-    else {
-      Msg::Error("%d-dimensional physical group %d does not exist",
-                 dim, basisPhysicals.at(i));
-      return 0;
-    }
+    curBasis.push_back(Chain<int>(m, basisPhysicals.at(i)));
   }
-  std::vector<std::string> curBasisNames2;
-  for(unsigned int i = 0; i < basisPhysicals2.size(); i++){
-    curBasisNames2.push_back(m->getPhysicalName(dim, basisPhysicals2.at(i)));
-    it = groups[dim].find(basisPhysicals2.at(i));
-    if(it != groups[dim].end()){
-      std::vector<GEntity*> physicalGroup = it->second;
-      for(unsigned int k = 0; k < physicalGroup.size(); k++){
-        basisEntities2.push_back(physicalGroup.at(k));
-        break;
-      }
-    }
-    else {
-      Msg::Error("%d-dimensional physical group %d does not exist",
-                 dim, basisPhysicals2.at(i));
-      return 0;
-    }
+  if(curBasis.empty()) {
+    Msg::Error("No operated chains given");
+    return 0;
   }
+  int dim = curBasis.at(0).getDim();
 
-  std::vector<Chain<int> > curBasis;
-  createChains(dim, m, basisEntities, curBasisNames, curBasis);
   std::vector<Chain<int> > curBasis2;
-  createChains(dim, m, basisEntities2, curBasisNames2, curBasis2);
-
+  for(unsigned int i = 0; i < basisPhysicals2.size(); i++){
+    curBasis2.push_back(Chain<int>(m, basisPhysicals2.at(i)));
+  }
 
   if(!curBasis2.empty()) {
     rows = curBasis2.size();
@@ -295,6 +239,9 @@ PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
   if(!curBasis2.empty()) {
     Msg::Info("Computing new basis %d-chains such that the incidence matrix of %d-chain bases %s and %s is the indentity matrix",
               dim, dim, opString1.c_str(), opString2.c_str());
+    int det = detIntegerMatrix(matrix);
+    if(det != 1 && det != -1)
+      Msg::Warning("Incidence matrix is not unimodular (det = %d)", det);
     if(!invertIntegerMatrix(matrix)) return 0;
     for(int i = 0; i < rows; i++)
       for(int j = 0; j < cols; j++)
@@ -305,7 +252,7 @@ PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
     if(rows == cols) {
       int det = detIntegerMatrix(matrix);
       if(det != 1 && det != -1)
-      Msg::Warning("Transformation matrix is not unimodular (det = %d)", det);
+        Msg::Warning("Transformation matrix is not unimodular (det = %d)", det);
     }
     for(int i = 0; i < rows; i++)
       for(int j = 0; j < cols; j++)
@@ -322,27 +269,8 @@ PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
   if(!tracePhysicals.empty()) {
     Msg::Info("Taking trace of result %d-chains to domain %s",
               dim, traceString.c_str());
-    std::vector<GEntity*> traceEntities;
-    for(unsigned int i = 0; i < tracePhysicals.size(); i++){
-      bool found = false;
-      for(int j = 0; j < 4; j++){
-        it = groups[j].find(tracePhysicals.at(i));
-        if(it != groups[j].end()){
-          found = true;
-          std::vector<GEntity*> physicalGroup = it->second;
-          for(unsigned int k = 0; k < physicalGroup.size(); k++){
-            traceEntities.push_back(physicalGroup.at(k));
-          }
-        }
-      }
-      if(!found) {
-        Msg::Error("Physical group %d does not exist",
-                   tracePhysicals.at(i));
-        return 0;
-      }
-    }
     for(unsigned int i = 0; i < newBasis.size(); i++)
-      newBasis.at(i) = newBasis.at(i).getTrace(traceEntities);
+      newBasis.at(i) = newBasis.at(i).getTrace(m, tracePhysicals);
     ElemChain::clearVertexCache();
   }
 
diff --git a/Plugin/HomologyPostProcessing.h b/Plugin/HomologyPostProcessing.h
index bcaa5e109ef815b4b8f66f9f9e2e89c476ff61f4..6c585196fcde0143620f5590d604990d3f4f4a42 100644
--- a/Plugin/HomologyPostProcessing.h
+++ b/Plugin/HomologyPostProcessing.h
@@ -36,10 +36,6 @@ class GMSH_HomologyPostProcessingPlugin : public GMSH_PostPlugin
   StringXString *getOptionStr(int iopt);
   PView *execute(PView *);
   bool parseStringOpt(int stringOpt, std::vector<int>& intList);
-  void createChains(int dim, GModel* m,
-                    const std::vector<GEntity*>& chainEntities,
-                    const std::vector<std::string>& chainNames,
-                    std::vector<Chain<int> >& chains);
   bool invertIntegerMatrix(std::vector<int>& matrix);
   int detIntegerMatrix(std::vector<int>& matrix);
 };