diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 8aa0c61eead299898d62e67cfb9590360dd35b31..dc45ef6c85f4515bb18c51f1b14b2889e81a06a0 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -96,6 +96,7 @@ void PrintUsage(const char *name)
   Msg::Direct("  -fontsize int         Specify the font size for the GUI");
   Msg::Direct("  -theme string         Specify FLTK GUI theme");
   Msg::Direct("  -display string       Specify display");
+  Msg::Direct("  -compound_only        Hide underlying surfaces/edges of compounds");
 #endif
   Msg::Direct("Other options:");      
   Msg::Direct("  -                     Parse input files, then exit");
@@ -733,6 +734,10 @@ void GetOptions(int argc, char *argv[])
         else
           Msg::Fatal("Missing argument");
       }
+      else if(!strcmp(argv[i] + 1, "compound_only")) {
+        CTX::instance()->compoundOnly = 1;
+        i++;
+      }
 #endif
 #if defined(__APPLE__)
       else if(!strncmp(argv[i] + 1, "psn", 3)) {
diff --git a/Common/Context.cpp b/Common/Context.cpp
index 1c4fcba6689c1df885e642fa622a6583e92cbe1c..36681e5500eb363a841d56d2cf7270f97b81b4d0 100644
--- a/Common/Context.cpp
+++ b/Common/Context.cpp
@@ -62,6 +62,7 @@ CTX::CTX()
 #endif
   forcedBBox = 0;
   hideUnselected = 0;
+  compoundOnly = 1;
   numWindows = numTiles = 1;
   deltaFontSize = 0;
   recentFiles.resize(5);
diff --git a/Common/Context.h b/Common/Context.h
index e21dfda744a41de395df8c76dd6c4a7b65344cce..f2c7acc17b311bd1ed40d7009cf633bda330327e 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -192,6 +192,8 @@ class CTX {
   int printing;
   // hide all unselected entities?
   int hideUnselected;
+  // hide underlying curves and surfaces of compounds (makes work a lot easier)
+  int compoundOnly;
   // temporary storage of rotation, translation, scale (until the GUI
   // is ready)
   double tmpRotation[3], tmpTranslation[3], tmpScale[3], tmpQuaternion[4];
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index d04e042aa1ec19cc66fef19765c81189ceb68a7a..a339b42a34bf84aec37bbf70f91d053febf47561 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -162,10 +162,19 @@ SOrientedBoundingBox GEdge::getOBB()
 
 void GEdge::setVisibility(char val, bool recursive)
 {
-  GEntity::setVisibility(val);
-  if(recursive){
-    if(v0) v0->setVisibility(val);
-    if(v1) v1->setVisibility(val);
+  if (getCompound() && CTX::instance()->compoundOnly) {
+    // use visibility info of compound edge if this edge belongs to it 
+    GEntity::setVisibility(0);
+    if(v0) v0->setVisibility(0);
+    if(v1) v1->setVisibility(0);
+    if(getCompound()->v0) getCompound()->v0->setVisibility(1);
+    if(getCompound()->v1) getCompound()->v1->setVisibility(1);
+  } else {
+    GEntity::setVisibility(val);
+    if(recursive){
+      if(v0) v0->setVisibility(val);
+      if(v1) v1->setVisibility(val);
+    }
   }
 }
 
diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index 2862238896f01278318d4b05325e7a4b98a4805c..a2c47fc78129880540f92704fd538e6b0c8cb91a 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -34,7 +34,7 @@ class GEdgeCompound : public GEdge {
   virtual double curvature(double t) const;
   virtual int minimumMeshSegments() const;
   virtual int minimumDrawSegments() const;
-  std::vector<GEdge*>  getEdgesOfCompound() const { return _compound; }
+  std::vector<GEdge*>  getCompounds() const { return _compound; }
 };
 
 #endif
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 7321425206845b7cfbad23d64e1b4c94b7b212aa..639273692639750885a37afd38a5a48b9229de82 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -242,12 +242,19 @@ std::list<GVertex*> GFace::vertices() const
 
 void GFace::setVisibility(char val, bool recursive)
 {
-  GEntity::setVisibility(val);
-  if(recursive){
-    std::list<GEdge*>::iterator it = l_edges.begin();
-    while (it != l_edges.end()){
-      (*it)->setVisibility(val, recursive);
-      ++it;
+  if (getCompound() && CTX::instance()->compoundOnly) {
+    GEntity::setVisibility(0);
+    for (std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); ++it)
+      (*it)->setVisibility(0, true);
+    std::list<GEdge*> edgesComp = getCompound()->edges();
+    //show edges of the compound surface
+    for (std::list<GEdge*>::iterator it = edgesComp.begin(); it != edgesComp.end(); ++it)
+      (*it)->setVisibility(1, true);
+  } else {
+    GEntity::setVisibility(val);
+    if(recursive){
+      for (std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); ++it)
+        (*it)->setVisibility(val, recursive);
     }
   }
 }
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 280e5554d85b8c491e682567d812436704da6def..55e516e611643c546655fbcafbd1bc8314d35aa4 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -984,7 +984,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
 
       //compute local parameter on Edge
       gec->getLocalParameter(tGlob,iEdge,tLoc);
-      std::vector<GEdge*> gev = gec->getEdgesOfCompound();
+      std::vector<GEdge*> gev = gec->getCompounds();
       GEdge *ge = gev[iEdge];
        
       //left and right vertex of the Edge
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 6005e3357e71a7f17c230f3d09882eea3148b478..5c5cc990dffa19ca095fa21500333df34f34546c 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -474,6 +474,23 @@ void GModel::setSelection(int val)
   }
 }
 
+void GModel::updateUpperTopology()
+{
+  for(eiter it = firstEdge(); it != lastEdge(); ++it) {
+    if(!(*it)->getCompound() && !edgesUpper.count((*it)))
+      edgesUpper.insert((*it));
+  }
+  for(fiter it = firstFace(); it != lastFace(); ++it) {
+    if(!(*it)->getCompound() && !facesUpper.count((*it)))
+      facesUpper.insert((*it));
+  }
+  for(riter it = firstRegion(); it != lastRegion(); ++it) {
+    if(!(*it)->getCompound() && !regionsUpper.count((*it)))
+      regionsUpper.insert((*it));
+  }
+}
+
+
 SBoundingBox3d GModel::bounds(bool aroundVisible)
 {
   std::vector<GEntity*> entities;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 44a2b2ffaceceaf018e96961332c89df6a6d5603..3137fb8f420e5ad588cd3e32b7df32024455ecf3 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -116,6 +116,12 @@ class GModel
   std::set<GFace*, GEntityLessThan> faces;
   std::set<GEdge*, GEntityLessThan> edges;
   std::set<GVertex*, GEntityLessThan> vertices;
+  
+  // represents uppermost topology level, when using compounds
+  // same as normal regions, edges, vertices but without entities contained in compounds
+  std::set<GRegion*, GEntityLessThan> regionsUpper;
+  std::set<GFace*, GEntityLessThan> facesUpper;
+  std::set<GEdge*, GEntityLessThan> edgesUpper;
 
   // map between the pair <dimension, elementary or physical number>
   // and an optional associated name
@@ -157,6 +163,9 @@ class GModel
   OCC_Internals *getOCCInternals(){ return _occ_internals; }
   FM_Internals *getFMInternals() { return _fm_internals; }
   ACIS_Internals *getACISInternals(){ return _acis_internals; }
+  
+  // if model has been loaded or changed fill facesUpper, edgesUpper
+  void updateUpperTopology();
 
   // access characteristic length (mesh size) fields
   FieldManager *getFields(){ return _fields; }
@@ -190,13 +199,13 @@ class GModel
   typedef std::set<GVertex*, GEntityLessThan>::iterator viter;
 
   // get an iterator initialized to the first/last entity in this model
-  riter firstRegion() { return regions.begin(); }
-  fiter firstFace() { return faces.begin(); }
-  eiter firstEdge() { return edges.begin(); }
+  riter firstRegion(const bool upper = false) { return upper ? regionsUpper.begin() : regions.begin(); }
+  fiter firstFace(const bool upper = false)   { return upper ? facesUpper.begin() : faces.begin(); }
+  eiter firstEdge(const bool upper = false)   { return upper ? edgesUpper.begin() : edges.begin(); }
   viter firstVertex() { return vertices.begin(); }
-  riter lastRegion() { return regions.end(); }
-  fiter lastFace() { return faces.end(); }
-  eiter lastEdge() { return edges.end(); }
+  riter lastRegion(const bool upper = false)  { return upper ? regionsUpper.end() : regions.end(); }
+  fiter lastFace(const bool upper = false)    { return upper ? facesUpper.end() : faces.end(); }
+  eiter lastEdge(const bool upper = false)    { return upper ? edgesUpper.end() : edges.end(); }
   viter lastVertex() { return vertices.end(); }
 
   // find the entity with the given tag
diff --git a/Geo/GModelIO_ACIS.cpp b/Geo/GModelIO_ACIS.cpp
index cd39c20e4372f73c54e5a437b0e01cf95e666cec..29d162bda16e5709b55a4d85263845728f553fe3 100644
--- a/Geo/GModelIO_ACIS.cpp
+++ b/Geo/GModelIO_ACIS.cpp
@@ -247,6 +247,7 @@ void ACIS_Internals::loadSAT(std::string fileName, GModel *gm)
       }
     }
   }
+  gm->updateUpperTopology();
 }
 
 int GModel::readACISSAT(const std::string &fn)
diff --git a/Geo/GModelIO_Fourier.cpp b/Geo/GModelIO_Fourier.cpp
index 2f63b4e1c762ac386a0db9ef2d3e95ee601f198b..7c44dc3b3b093bda0da1dff6211fbbc03250c7e7 100644
--- a/Geo/GModelIO_Fourier.cpp
+++ b/Geo/GModelIO_Fourier.cpp
@@ -97,6 +97,7 @@ void FM_Internals::buildGModel(FM::Reader* reader, GModel* model)
 {
   for (int i = 0; i < reader->GetNumPatches(); i++)
     makeGFace(reader->GetPatch(i), model);
+  model->updateUpperTopology();
 }
 
 void GModel::_createFMInternals()
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 6afe76ae17fc009cef327d4d2ecb83a7bb09e59b..fcb320998410065f58aff8b4927efee88e2c7fa1 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -177,50 +177,53 @@ int GModel::importGEOInternals()
         std::list<GFace*> comp;
         for(unsigned int j = 0; j < s->compound.size(); j++){
           GFace *gf = getFaceByTag(s->compound[j]);
-          if(gf) comp.push_back(gf);
+          if(gf)
+            comp.push_back(gf);
         }
-	std::list<GEdge*> U0;
+        std::list<GEdge*> U0;
         for(int j = 0; j < 4; j++){
           for(unsigned int k = 0; k < s->compoundBoundary[j].size(); k++){
             GEdge *ge = getEdgeByTag(s->compoundBoundary[j][k]);
             if(ge) U0.push_back(ge);
           }
         }
- 	int param = CTX::instance()->mesh.remeshParam;
-	GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
-	if (param == 1) typ =  GFaceCompound::CONFORMAL;
-	if (param == 2) typ =  GFaceCompound::RBF;
-	int algo = CTX::instance()->mesh.remeshAlgo;
-        f = new GFaceCompound(this, s->Num, comp, U0, typ, algo);
+        int param = CTX::instance()->mesh.remeshParam;
+        GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
+        if (param == 1) typ =  GFaceCompound::CONFORMAL;
+        if (param == 2) typ =  GFaceCompound::RBF;
+        int algo = CTX::instance()->mesh.remeshAlgo;
+              f = new GFaceCompound(this, s->Num, comp, U0, typ, algo);
 
-	f->meshAttributes.recombine = s->Recombine;
-	f->meshAttributes.recombineAngle = s->RecombineAngle;
-	f->meshAttributes.Method = s->Method;
-	f->meshAttributes.extrude = s->Extrude;
+        f->meshAttributes.recombine = s->Recombine;
+        f->meshAttributes.recombineAngle = s->RecombineAngle;
+        f->meshAttributes.Method = s->Method;
+        f->meshAttributes.extrude = s->Extrude;
         add(f);
-
-	if(s->EmbeddedCurves){
-	  for(int i = 0; i < List_Nbr(s->EmbeddedCurves); i++){
-	    Curve *c;
-	    List_Read(s->EmbeddedCurves, i, &c);
-	    GEdge *e = getEdgeByTag(abs(c->Num));
-	    if(e)
-	      f->addEmbeddedEdge(e);
-	    else
-	      Msg::Error("Unknown curve %d", c->Num);
-	  }
-	}
-	if(s->EmbeddedPoints){
-	  for(int i = 0; i < List_Nbr(s->EmbeddedPoints); i++){
-	    Vertex *v;
-	    List_Read(s->EmbeddedPoints, i, &v);
-	    GVertex *gv = getVertexByTag(v->Num);
-	    if(gv)
-	      f->addEmbeddedVertex(gv);
-	    else
-	      Msg::Error("Unknown point %d", v->Num);
-	  }
-	}
+        if (CTX::instance()->compoundOnly)
+          for (std::list<GFace*>::iterator it = comp.begin(); it != comp.end(); ++it)
+            (*it)->setVisibility(0,true);
+        if(s->EmbeddedCurves){
+          for(int i = 0; i < List_Nbr(s->EmbeddedCurves); i++){
+            Curve *c;
+            List_Read(s->EmbeddedCurves, i, &c);
+            GEdge *e = getEdgeByTag(abs(c->Num));
+            if(e)
+              f->addEmbeddedEdge(e);
+            else
+              Msg::Error("Unknown curve %d", c->Num);
+          }
+        }
+        if(s->EmbeddedPoints){
+          for(int i = 0; i < List_Nbr(s->EmbeddedPoints); i++){
+            Vertex *v;
+            List_Read(s->EmbeddedPoints, i, &v);
+            GVertex *gv = getVertexByTag(v->Num);
+            if(gv)
+              f->addEmbeddedVertex(gv);
+            else
+              Msg::Error("Unknown point %d", v->Num);
+          }
+        }
       }
       else if(!f){
         f = new gmshFace(this, s);
@@ -279,6 +282,8 @@ int GModel::importGEOInternals()
     }
   }
 
+  updateUpperTopology();
+
   Msg::Debug("Gmsh model (GModel) imported:");
   Msg::Debug("%d Vertices", vertices.size());
   Msg::Debug("%d Edges", edges.size());
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 4d80c00183f796e8062a6229f80778ed688d729b..b7ae3887d2eb64e68ba634f6814b446a2522824b 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -707,6 +707,8 @@ void OCC_Internals::buildGModel(GModel *model)
     if (!getOCCRegionByNativePtr(model, TopoDS::Solid(somap(i))))
       model->add(new OCCRegion(model, TopoDS::Solid(somap(i)), num));
   }
+  
+  model->updateUpperTopology();
 }
 
 void addSimpleShapes(TopoDS_Shape theShape, TopTools_ListOfShape &theList)
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index 77acd62130f7ccb2988fb03368826259b34bad02..0770c3007bbfc98711aa300d166a854e4505e4f5 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -19,13 +19,15 @@ class MPrism;
 class MPyramid;
 class MPolyhedron;
 class ExtrudeParams;
+class GRegionCompound;
 
 // A model region.
 class GRegion : public GEntity {
  protected:
   std::list<GFace*> l_faces;
   std::list<int> l_dirs;
-
+  GRegionCompound *compound; // this model ede belongs to a compound 
+  
   // replace faces (for gluing) for specific modelers, we have to
   // re-create internal data ...
   virtual void replaceFacesInternal (std::list<GFace*> &) {}
@@ -89,6 +91,10 @@ class GRegion : public GEntity {
 
   // reset the mesh attributes to default values
   virtual void resetMeshAttributes();
+  
+  // compound
+  void setCompound(GRegionCompound *grc) { compound = grc; }
+  GRegionCompound *getCompound() const { return compound; }
 
   struct {
     // is this surface meshed using a transfinite interpolation
diff --git a/Geo/GRegionCompound.cpp b/Geo/GRegionCompound.cpp
index a89a37a1081b4401333ce69d78acfd857aca0391..97218114966d5c510b869d349073398a1832b3f4 100644
--- a/Geo/GRegionCompound.cpp
+++ b/Geo/GRegionCompound.cpp
@@ -15,6 +15,15 @@
 GRegionCompound::GRegionCompound(GModel *m, int tag, std::vector<GRegion*> &compound)
   : GRegion(m, tag), _compound(compound)
 {
+  
+  for (unsigned int i = 0; i < _compound.size(); i++){
+    if(!_compound[i]){
+      Msg::Error("Incorrect region in compound region %d\n", tag);
+      Msg::Exit(1);
+    }
+  }
+  for (unsigned int i = 0; i < _compound.size(); i++)
+    _compound[i]->setCompound(this);
   getBoundingFaces();
 }
 
diff --git a/Graphics/drawGeom.cpp b/Graphics/drawGeom.cpp
index f4634cffe6f078e559012be70ce151274e16fb43..3a8625e8a937b1d8eee3c933bb8e65ac1e88d643 100644
--- a/Graphics/drawGeom.cpp
+++ b/Graphics/drawGeom.cpp
@@ -97,6 +97,7 @@ class drawGEdge {
     if(e->geomType() == GEntity::DiscreteCurve) return;
     if(e->geomType() == GEntity::PartitionCurve) return;
     if(e->geomType() == GEntity::BoundaryLayerCurve) return;
+    if(e->geomType() == GEntity::CompoundCurve) return; 
     
     bool select = (_ctx->render_mode == drawContext::GMSH_SELECT && 
                    e->model() == GModel::current());
@@ -236,8 +237,6 @@ class drawGFace {
   }
   void _drawParametricGFace(GFace *f)
   {
-    if (f->geomType() == GEntity::CompoundSurface) return;
-
     if(CTX::instance()->geom.surfaceType > 0)
       f->fillVertexArray(f->geomType() == GEntity::ProjectionFace);
 
@@ -248,9 +247,12 @@ class drawGFace {
 
     if(CTX::instance()->geom.surfaces){
       if(CTX::instance()->geom.surfaceType > 0 && f->va_geom_triangles){
+        bool selected = false;
+        if (f->getSelection() || (f->getCompound() && f->getCompound()->getSelection()))
+          selected = true;
         _drawVertexArray
           (f->va_geom_triangles, CTX::instance()->geom.light, 
-           (f->geomType() == GEntity::ProjectionFace) ? true : f->getSelection(), 
+           (f->geomType() == GEntity::ProjectionFace) ? true : selected, 
            (f->geomType() == GEntity::ProjectionFace) ? 
            CTX::instance()->color.geom.projection : 
            CTX::instance()->color.geom.selection);
@@ -319,6 +321,9 @@ class drawGFace {
       f->buildRepresentationCross();
 
     if(CTX::instance()->geom.surfaces) {
+      bool selected = false;
+      if (f->getSelection() || (f->getCompound() && f->getCompound()->getSelection()))
+          selected = true;
       if(CTX::instance()->geom.surfaceType > 0 && f->va_geom_triangles){
         _drawVertexArray(f->va_geom_triangles, CTX::instance()->geom.light, 
                          f->getSelection(), CTX::instance()->color.geom.selection);
@@ -371,14 +376,32 @@ class drawGFace {
                        CTX::instance()->geom.light);
     }
   }
+  
+  void _drawCompoundGFace(GFace *f, bool visible = false, bool selected = false)
+  {
+    GFaceCompound *fc = (GFaceCompound*) f;
+    std::list<GFace*> faces = fc->getCompounds();
+    for (std::list<GFace*>::iterator it = faces.begin(); it!=faces.end(); it++) {
+      if ((*it)->geomType() == GEntity::DiscreteSurface) continue;
+      if ((*it)->geomType() == GEntity::PartitionSurface) continue;
+      if ((*it)->geomType() == GEntity::BoundaryLayerSurface) continue;
+        
+      if((*it)->geomType() == GEntity::CompoundSurface)
+        _drawCompoundGFace((*it));
+      else if ((*it)->geomType() == GEntity::Plane)
+        _drawPlaneGFace((*it));
+      else
+        _drawParametricGFace((*it));
+    }
+  }
  
  public :
   drawGFace(drawContext *ctx) : _ctx(ctx) {}
   void operator () (GFace *f)
   {
     if(!f->getVisibility()) return;
-     if(f->geomType() == GEntity::DiscreteSurface) return;
-     if(f->geomType() == GEntity::PartitionSurface) return;
+    if(f->geomType() == GEntity::DiscreteSurface) return;
+    if(f->geomType() == GEntity::PartitionSurface) return;
     if(f->geomType() == GEntity::BoundaryLayerSurface) return;
 
     bool select = (_ctx->render_mode == drawContext::GMSH_SELECT && 
@@ -401,7 +424,9 @@ class drawGFace {
       glColor4ubv((GLubyte *) & CTX::instance()->color.geom.surface);
     }
 
-    if(f->geomType() == GEntity::Plane)
+    if(f->geomType() == GEntity::CompoundSurface)
+      _drawCompoundGFace(f);
+    else if(f->geomType() == GEntity::Plane)
       _drawPlaneGFace(f);
     else
       _drawParametricGFace(f);
@@ -484,12 +509,14 @@ void drawContext::drawGeom()
         std::for_each(m->firstVertex(), m->lastVertex(), drawGVertex(this));
       if(CTX::instance()->geom.lines || CTX::instance()->geom.linesNum || 
          CTX::instance()->geom.tangents)
-        std::for_each(m->firstEdge(), m->lastEdge(), drawGEdge(this));
+        std::for_each(m->firstEdge(true), m->lastEdge(true), drawGEdge(this));
       if(CTX::instance()->geom.surfaces || CTX::instance()->geom.surfacesNum ||
-         CTX::instance()->geom.normals)
-        std::for_each(m->firstFace(), m->lastFace(), drawGFace(this));
+         CTX::instance()->geom.normals) {
+        GModel::fiter test = m->firstFace(true); 
+        std::for_each(m->firstFace(true), m->lastFace(true), drawGFace(this));
+      }
       if(CTX::instance()->geom.volumes || CTX::instance()->geom.volumesNum)
-        std::for_each(m->firstRegion(), m->lastRegion(), drawGRegion(this));
+        std::for_each(m->firstRegion(true), m->lastRegion(true), drawGRegion(this));
     }
   }