diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 59117e31ac21ed95b9db5bf8384ded18ee85d08b..1d7865613f30c120c8582c4c10d1837ed2b04848 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -288,14 +288,13 @@ GVertex *GModel::getVertexByTag(int n) const
 std::vector<int> GModel::getEdgesByStringTag(const std::string tag)
 {
   std::vector<int> nums;
-  std::map<int, std::vector<GEntity*> > physicalGroups[4];
-  getPhysicalGroups(physicalGroups);
-  std::vector<GEntity*> allEdges = physicalGroups[1][this->getPhysicalNumber(1,tag)];
-  for ( std::vector<GEntity*>::iterator it = allEdges.begin(); it != allEdges.end(); it++){
+  std::map<int, std::vector<GEntity*> > physicalGroups;
+  getPhysicalGroups(1, physicalGroups);
+  std::vector<GEntity*> allEdges = physicalGroups[this->getPhysicalNumber(1, tag)];
+  for (std::vector<GEntity*>::iterator it = allEdges.begin(); it != allEdges.end(); it++){
     GEntity *ge = *it;
     nums.push_back(ge->tag());
   }
-
   return nums;
 }
 
@@ -403,7 +402,7 @@ bool GModel::noPhysicalGroups()
   return true;
 }
 
-void GModel::getPhysicalGroups(std::map<int, std::vector<GEntity*> > groups[4])
+void GModel::getPhysicalGroups(std::map<int, std::vector<GEntity*> > groups[4]) const
 {
   std::vector<GEntity*> entities;
   getEntities(entities);
@@ -419,6 +418,21 @@ void GModel::getPhysicalGroups(std::map<int, std::vector<GEntity*> > groups[4])
   }
 }
 
+void GModel::getPhysicalGroups(int dim, std::map<int, std::vector<GEntity*> > &groups) const
+{
+  std::vector<GEntity*> entities;
+  getEntities(entities, dim);
+  for(unsigned int i = 0; i < entities.size(); i++){
+    for(unsigned int j = 0; j < entities[i]->physicals.size(); j++){
+      // physicals can be stored with negative signs when the entity
+      // should be "reversed"
+      int p = std::abs(entities[i]->physicals[j]);
+      if(std::find(groups[p].begin(), groups[p].end(), entities[i]) == groups[p].end())
+        groups[p].push_back(entities[i]);
+    }
+  }
+}
+
 void GModel::deletePhysicalGroups()
 {
   std::vector<GEntity*> entities;
@@ -430,15 +444,13 @@ void GModel::deletePhysicalGroups()
 void GModel::deletePhysicalGroup(int dim, int num)
 {
   std::vector<GEntity*> entities;
-  getEntities(entities);
+  getEntities(entities, dim);
   for(unsigned int i = 0; i < entities.size(); i++){
-    if(dim == entities[i]->dim()){
-      std::vector<int> p;
-      for(unsigned int j = 0; j < entities[i]->physicals.size(); j++)
-        if(entities[i]->physicals[j] != num)
-          p.push_back(entities[i]->physicals[j]);
-      entities[i]->physicals = p;
-    }
+    std::vector<int> p;
+    for(unsigned int j = 0; j < entities[i]->physicals.size(); j++)
+      if(entities[i]->physicals[j] != num)
+        p.push_back(entities[i]->physicals[j]);
+    entities[i]->physicals = p;
   }
 }
 
@@ -470,14 +482,6 @@ int GModel::setPhysicalName(std::string name, int dim, int number)
 
 std::string GModel::getPhysicalName(int dim, int number) const
 {
-  //Emi debug here
-  // printf("getPhysName size %d \n", physicalNames.size());
-  // std::map<std::pair<int, int>, std::string>::iterator itt = physicalNames.begin();
-  // for (; itt != physicalNames.end(); itt++){
-  //   printf("name %s \n", itt->second.c_str());
-  //   printf("par (%d,%d) \n", itt->first.first, itt->first.second);
-  // }
-
   std::map<std::pair<int, int>, std::string>::const_iterator it =
     physicalNames.find(std::pair<int, int>(dim, number));
   if(it != physicalNames.end()) return it->second;
@@ -836,10 +840,10 @@ MVertex *GModel::getMeshVertexByTag(int n)
 void GModel::getMeshVerticesForPhysicalGroup(int dim, int num, std::vector<MVertex*> &v)
 {
   v.clear();
-  std::map<int, std::vector<GEntity*> > groups[4];
-  getPhysicalGroups(groups);
-  std::map<int, std::vector<GEntity*> >::const_iterator it = groups[dim].find(num);
-  if(it == groups[dim].end()) return;
+  std::map<int, std::vector<GEntity*> > groups;
+  getPhysicalGroups(dim, groups);
+  std::map<int, std::vector<GEntity*> >::const_iterator it = groups.find(num);
+  if(it == groups.end()) return;
   const std::vector<GEntity *> &entities = it->second;
   std::set<MVertex*> sv;
   for(unsigned int i = 0; i < entities.size(); i++){
@@ -1033,8 +1037,8 @@ void GModel::deleteMeshPartitions()
 }
 
 void GModel::store(std::vector<MVertex*> &vertices, int dim,
-          std::map<int, std::vector<MElement*> > &entityMap,
-          std::map<int, std::map<int, std::string> > &physicalMap)
+                   std::map<int, std::vector<MElement*> > &entityMap,
+                   std::map<int, std::map<int, std::string> > &physicalMap)
 {
   _storeVerticesInEntities(vertices);
   _storeElementsInEntities(entityMap);
diff --git a/Geo/GModel.h b/Geo/GModel.h
index fc598af855e9563ac344ff851febf9c742200cd6..5b17c53ad6d01bd7c329f57b15279688e4b8cf97 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -292,7 +292,8 @@ class GModel
   bool noPhysicalGroups();
 
   // return all physical groups (one map per dimension: 0-D to 3-D)
-  void getPhysicalGroups(std::map<int, std::vector<GEntity*> > groups[4]);
+  void getPhysicalGroups(std::map<int, std::vector<GEntity*> > groups[4]) const;
+  void getPhysicalGroups(int dim, std::map<int, std::vector<GEntity*> > &groups) const;
 
   // delete physical groups in the model
   void deletePhysicalGroups();
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index f9dec983450f38da10ef1a731b2e8bea6d8dc733..10480f27339222307927b6c4e8e63328503fd8e7 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -20,7 +20,6 @@
 #include "Geo.h"
 #include "GmshDefines.h"
 
-
 GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
 {
   int num =  gm->getMaxElementaryNumber(0) + 1;
@@ -164,7 +163,6 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
 
 GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > faces)
 {
-
   //create surface loop
   int nLoops = faces.size();
   std::vector<SurfaceLoop *> vecLoops;
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 636f88493987bc6c3b01a8852f7201f4dfb8b28e..f5751aa0a542bc1434f91d5f774d9da956bd10e5 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -1087,9 +1087,9 @@ gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv)
 gLevelsetYarn::gLevelsetYarn(int dim, int phys, double minA, double majA, int type, int tag)
   : gLevelsetPrimitive(tag), minorAxis(minA), majorAxis(majA), typeLs(type)
 {
-  std::map<int, std::vector<GEntity*> > groups[4];
-  GModel::current()->getPhysicalGroups(groups);
-  entities = groups[dim][phys];
+  std::map<int, std::vector<GEntity*> > groups;
+  GModel::current()->getPhysicalGroups(dim, groups);
+  entities = groups[phys];
   if(!entities.size())
     printf("No physical %d found for levelset yarn!\n", phys);
 }
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 12ed003a053f47547513305ab0c1568e33e40bb9..d38626a4a1bc6c823833ce719e56eac6c407c0c7 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -10490,10 +10490,10 @@ yyreduce:
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[0].find((int)num);
-          if(it != groups[0].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(0, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add((yyval.l), &d);
@@ -10523,10 +10523,10 @@ yyreduce:
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[1].find((int)num);
-          if(it != groups[1].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(1, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add((yyval.l), &d);
@@ -10556,10 +10556,10 @@ yyreduce:
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[2].find((int)num);
-          if(it != groups[2].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(2, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add((yyval.l), &d);
@@ -10589,10 +10589,10 @@ yyreduce:
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[3].find((int)num);
-          if(it != groups[3].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(3, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add((yyval.l), &d);
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index a8bcd7457d1ee624e74931a31bc61b2f518e195d..6892d77115583d354e6d3af8f00dd1446b22f330 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -4813,10 +4813,10 @@ FExpr_Multi :
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[0].find((int)num);
-          if(it != groups[0].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(0, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add($$, &d);
@@ -4842,10 +4842,10 @@ FExpr_Multi :
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[1].find((int)num);
-          if(it != groups[1].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(1, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add($$, &d);
@@ -4871,10 +4871,10 @@ FExpr_Multi :
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[2].find((int)num);
-          if(it != groups[2].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(2, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add($$, &d);
@@ -4900,10 +4900,10 @@ FExpr_Multi :
           }
         }
         else{
-          std::map<int, std::vector<GEntity*> > groups[4];
-          GModel::current()->getPhysicalGroups(groups);
-          std::map<int, std::vector<GEntity*> >::iterator it = groups[3].find((int)num);
-          if(it != groups[3].end()){
+          std::map<int, std::vector<GEntity*> > groups;
+          GModel::current()->getPhysicalGroups(3, groups);
+          std::map<int, std::vector<GEntity*> >::iterator it = groups.find((int)num);
+          if(it != groups.end()){
             for(unsigned j = 0; j < it->second.size(); j++){
               double d = it->second[j]->tag();
               List_Add($$, &d);
diff --git a/Plugin/MeshSubEntities.cpp b/Plugin/MeshSubEntities.cpp
index 52b84bef71d05b9ca0ca3b8367bd5127f0336ce4..dd67f4bc698a45bee9e921a6203af9f7b79913b5 100644
--- a/Plugin/MeshSubEntities.cpp
+++ b/Plugin/MeshSubEntities.cpp
@@ -65,9 +65,9 @@ PView *GMSH_MeshSubEntitiesPlugin::execute(PView *view)
   }
 
   GModel *m = GModel::current();
-  std::map<int, std::vector<GEntity*> > groups[4];
-  m->getPhysicalGroups(groups);
-  std::vector<GEntity*> entities = groups[inputdim][inputphysical];
+  std::map<int, std::vector<GEntity*> > groups;
+  m->getPhysicalGroups(inputdim, groups);
+  std::vector<GEntity*> entities = groups[inputphysical];
 
   if(entities.empty()){
     Msg::Error("Physical group %d (dimension %d) is empty", inputphysical, inputdim);