diff --git a/Common/Context.h b/Common/Context.h
index d218da782dfa5763a5f91f334eac84a61fd1906b..c2a049ba5adf5ab6497f182ce3217cfc52558e17 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -226,7 +226,8 @@ class CTX {
     int fileFormat, epsQuality, epsCompress, epsPS3Shading;
     int epsOcclusionCulling, epsBestRoot;
     double epsLineWidthFactor, epsPointSizeFactor;
-    int jpegQuality, jpegSmoothing, geoLabels, text, texAsEquation;
+    int jpegQuality, jpegSmoothing, geoLabels, geoOnlyPhysicals;
+    int text, texAsEquation;
     int gifDither, gifSort, gifInterlace, gifTransparent;
     int posElementary, posElement, posGamma, posEta, posRho, posDisto;
     int compositeWindows, deleteTmpFiles, background;
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index bb2d0eaa3bfe022145302b69f88122b2c3447bfe..72a203a99fd5d67440b7319dfeb28405230e8084 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -290,7 +290,8 @@ void CreateOutputFile(std::string fileName, int format, bool redraw)
     break;
 
   case FORMAT_GEO:
-    GModel::current()->writeGEO(fileName, CTX::instance()->print.geoLabels);
+    GModel::current()->writeGEO(fileName, CTX::instance()->print.geoLabels,
+                                CTX::instance()->print.geoOnlyPhysicals);
     break;
 
   case FORMAT_BREP:
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 45565506b8b9f3a7b9345b966e9a3960b845e5dc..7f8027f8464ab2ea4caeb2969811a9a547b8fdbb 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1416,6 +1416,8 @@ StringXNumber PrintOptions_Number[] = {
 
   { F|O, "GeoLabels" , opt_print_geo_labels , 1. ,
     "Save labels in unrolled Gmsh geometries" },
+  { F|O, "GeoOnlyPhysicals" , opt_print_geo_only_physicals , 1. ,
+    "Only save entities that belong to physical groups" },
   { F|O, "GifDither" , opt_print_gif_dither , 0. ,
     "Apply dithering to GIF output" },
   { F|O, "GifInterlace" , opt_print_gif_interlace , 0. ,
diff --git a/Common/GmshSocket.h b/Common/GmshSocket.h
index 22673e2035c9cc9f1a9e2aeb58ddb34f0c3456a8..ee7671f8f624002746c9189ad3bc793b2e2f7fcb 100644
--- a/Common/GmshSocket.h
+++ b/Common/GmshSocket.h
@@ -320,7 +320,8 @@ class GmshServer : public GmshSocket{
   virtual ~GmshServer(){}
   virtual int NonBlockingSystemCall(const char *str) = 0;
   virtual int NonBlockingWait(int socket, double waitint, double timeout) = 0;
-
+  // start the client by launching "command" (command is supposed to contain
+  // '%s' where the socket name should appear)
   int Start(const char *command, const char *sockname, double timeout)
   {
     if(!sockname) throw "Invalid (null) socket name";
@@ -386,9 +387,9 @@ class GmshServer : public GmshSocket{
       }
     }
 
-    char cmd[1024];
     if(command && strlen(command)){
-      sprintf(cmd,command,_sockname.c_str());
+      char cmd[1024];
+      sprintf(cmd, command, _sockname.c_str());
       NonBlockingSystemCall(cmd); // starts the solver
     }
     else{
diff --git a/Common/Options.cpp b/Common/Options.cpp
index a300ef8b411c76a483676c7943cac4de2efe20dc..c9a6280a841c07a091c5075482f80a8c6246a25f 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -7874,6 +7874,13 @@ double opt_print_geo_labels(OPT_ARGS_NUM)
   return CTX::instance()->print.geoLabels;
 }
 
+double opt_print_geo_only_physicals(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX::instance()->print.geoOnlyPhysicals = (int)val;
+  return CTX::instance()->print.geoOnlyPhysicals;
+}
+
 double opt_print_pos_elementary(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 7a792231c583c6033fcfd314026067bdf1780e15..00597f506455e774a68b520ce1da46e2ba016e4c 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -592,6 +592,7 @@ double opt_print_eps_point_size_factor(OPT_ARGS_NUM);
 double opt_print_jpeg_quality(OPT_ARGS_NUM);
 double opt_print_jpeg_smoothing(OPT_ARGS_NUM);
 double opt_print_geo_labels(OPT_ARGS_NUM);
+double opt_print_geo_only_physicals(OPT_ARGS_NUM);
 double opt_print_pos_elementary(OPT_ARGS_NUM);
 double opt_print_pos_element(OPT_ARGS_NUM);
 double opt_print_pos_gamma(OPT_ARGS_NUM);
diff --git a/Fltk/fileDialogs.cpp b/Fltk/fileDialogs.cpp
index f4b145f79392bef71c05fb4c4b372ea0fddc06e3..f288d7d1bf0735aaef4e4cc58e17d0df0f068527 100644
--- a/Fltk/fileDialogs.cpp
+++ b/Fltk/fileDialogs.cpp
@@ -722,27 +722,31 @@ int geoFileDialog(const char *name)
 {
   struct _geoFileDialog{
     Fl_Window *window;
-    Fl_Check_Button *b;
+    Fl_Check_Button *b[2];
     Fl_Button *ok, *cancel;
   };
   static _geoFileDialog *dialog = NULL;
 
   if(!dialog){
     dialog = new _geoFileDialog;
-    int h = 3 * WB + 2 * BH, w = 2 * BB + 3 * WB, y = WB;
+    int h = 3 * WB + 3 * BH, w = 2 * BB + 3 * WB, y = WB;
     dialog->window = new Fl_Double_Window(w, h, "GEO Options");
     dialog->window->box(GMSH_WINDOW_BOX);
     dialog->window->set_modal();
-    dialog->b = new Fl_Check_Button
+    dialog->b[0] = new Fl_Check_Button
       (WB, y, 2 * BB + WB, BH, "Save physical group labels"); y += BH;
-    dialog->b->type(FL_TOGGLE_BUTTON);
+    dialog->b[0]->type(FL_TOGGLE_BUTTON);
+    dialog->b[1] = new Fl_Check_Button
+      (WB, y, 2 * BB + WB, BH, "Only save physical entities"); y += BH;
+    dialog->b[1]->type(FL_TOGGLE_BUTTON);
     dialog->ok = new Fl_Return_Button(WB, y + WB, BB, BH, "OK");
     dialog->cancel = new Fl_Button(2 * WB + BB, y + WB, BB, BH, "Cancel");
     dialog->window->end();
     dialog->window->hotspot(dialog->window);
   }
 
-  dialog->b->value(CTX::instance()->print.geoLabels ? 1 : 0);
+  dialog->b[0]->value(CTX::instance()->print.geoLabels ? 1 : 0);
+  dialog->b[1]->value(CTX::instance()->print.geoOnlyPhysicals ? 1 : 0);
   dialog->window->show();
 
   while(dialog->window->shown()){
@@ -751,7 +755,8 @@ int geoFileDialog(const char *name)
       Fl_Widget* o = Fl::readqueue();
       if (!o) break;
       if (o == dialog->ok) {
-        opt_print_geo_labels(0, GMSH_SET | GMSH_GUI, dialog->b->value() ? 1 : 0);
+        opt_print_geo_labels(0, GMSH_SET | GMSH_GUI, dialog->b[0]->value() ? 1 : 0);
+        opt_print_geo_only_physicals(0, GMSH_SET | GMSH_GUI, dialog->b[1]->value() ? 1 : 0);
         CreateOutputFile(name, FORMAT_GEO);
         dialog->window->hide();
         return 1;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 257f6da1f110f03ce9fe8e1386095571fe4ae48a..18d062f313ad2c41aa820050fb3e20def52a8ebe 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -509,7 +509,8 @@ class GModel
   // Gmsh native CAD format (readGEO is static, since it can create
   // multiple models)
   static int readGEO(const std::string &name);
-  int writeGEO(const std::string &name, bool printLabels=true);
+  int writeGEO(const std::string &name, bool printLabels=true,
+               bool onlyPhysicals=false);
   int importGEOInternals();
   int exportDiscreteGEOInternals();
 
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index e4df0f022a4ea39091b94d6b4725e5acbb961b8c..6ffa531e8252d1df9f49e144d1fd217e4ba1dcb3 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -374,7 +374,43 @@ class writePhysicalGroupGEO {
   }
 };
 
-int GModel::writeGEO(const std::string &name, bool printLabels)
+static bool skipRegion(GRegion *gr)
+{
+  if(gr->physicals.size()) return false;
+  return true;
+}
+
+static bool skipFace(GFace *gf)
+{
+  if(gf->physicals.size()) return false;
+  std::list<GRegion*> regions(gf->regions());
+  for(std::list<GRegion*>::iterator itr = regions.begin(); itr != regions.end(); itr++){
+    if(!skipRegion(*itr)) return false;
+  }
+  return true;
+}
+
+static bool skipEdge(GEdge *ge)
+{
+  if(ge->physicals.size()) return false;
+  std::list<GFace*> faces(ge->faces());
+  for(std::list<GFace*>::iterator itf = faces.begin(); itf != faces.end(); itf++){
+    if(!skipFace(*itf)) return false;
+  }
+  return true;
+}
+
+static bool skipVertex(GVertex *gv)
+{
+  if(gv->physicals.size()) return false;
+  std::list<GEdge*> edges(gv->edges());
+  for(std::list<GEdge*>::iterator ite = edges.begin(); ite != edges.end(); ite++){
+    if(!skipEdge(*ite)) return false;
+  }
+  return true;
+}
+
+int GModel::writeGEO(const std::string &name, bool printLabels, bool onlyPhysicals)
 {
   FILE *fp = fopen(name.c_str(), "w");
 
@@ -389,16 +425,24 @@ int GModel::writeGEO(const std::string &name, bool printLabels)
       meshSizeParameters.insert(std::make_pair(val, paramName.str()));
     }
   }
+
   for(viter it = firstVertex(); it != lastVertex(); it++){
     double val = (*it)->prescribedMeshSizeAtVertex();
-    (*it)->writeGEO(fp, meshSizeParameters[val]);
+    if(!onlyPhysicals || !skipVertex(*it))
+      (*it)->writeGEO(fp, meshSizeParameters[val]);
+  }
+  for(eiter it = firstEdge(); it != lastEdge(); it++){
+    if(!onlyPhysicals || !skipEdge(*it))
+      (*it)->writeGEO(fp);
+  }
+  for(fiter it = firstFace(); it != lastFace(); it++){
+    if(!onlyPhysicals || !skipFace(*it))
+      (*it)->writeGEO(fp);
+  }
+  for(riter it = firstRegion(); it != lastRegion(); it++){
+    if(!onlyPhysicals || !skipRegion(*it))
+      (*it)->writeGEO(fp);
   }
-  for(eiter it = firstEdge(); it != lastEdge(); it++)
-    (*it)->writeGEO(fp);
-  for(fiter it = firstFace(); it != lastFace(); it++)
-    (*it)->writeGEO(fp);
-  for(riter it = firstRegion(); it != lastRegion(); it++)
-    (*it)->writeGEO(fp);
 
   std::map<int, std::string> labels;
 #if defined(HAVE_PARSER)
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 004bcf76e3710a1ece054793597b688ef891f661..6729ff742dde1145b57896d47586ea30c042ab84 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -136,7 +136,7 @@ static MElement *createElementMSH(GModel *m, int num, int typeMSH, int physical,
     Msg::Error("Unknown type of element %d", typeMSH);
     return NULL;
   }
-  
+
 #if (FAST_ELEMENTS!=1)
   switch(e->getType()){
   case TYPE_PNT :
@@ -403,13 +403,13 @@ int GModel::readMSH(const std::string &name)
 
       if(!fgets(str, sizeof(str), fp)) return 0;
       int numElements;
-      
-      std::map<int, MElement*> elems;      
+
+      std::map<int, MElement*> elems;
       std::map<int, MElement*> parents;
       std::map<int, MElement*> domains;
       std::map<int, int> elemreg;
       std::map<int, int> elemphy;
-      
+
       std::set<MElement*> parentsOwned;
       sscanf(str, "%d", &numElements);
       Msg::Info("%d elements", numElements);
@@ -463,7 +463,7 @@ int GModel::readMSH(const std::string &name)
           }
           MElement *p = NULL;
           bool own = false;
-	  
+
 	  // search parent element
 	  if (parent)
 	  {
@@ -502,7 +502,7 @@ int GModel::readMSH(const std::string &name)
 	      }
 #endif
 	  }
-          
+
           // search domains
           MElement *doms[2] = {NULL, NULL};
 	  if(dom1)
@@ -517,7 +517,7 @@ int GModel::readMSH(const std::string &name)
 	    ite = elems.find(dom2);
 	    if (ite != elems.end())
 		doms[1] = ite->second;
-	    
+
             if(!doms[0]) Msg::Error("Domain element %d not found for element %d", dom1, num);
             if(dom2 && !doms[1]) Msg::Error("Domain element %d not found for element %d", dom2, num);
 #else
@@ -526,31 +526,31 @@ int GModel::readMSH(const std::string &name)
             if(dom2 && !doms[1]) Msg::Error("Domain element %d not found for element %d", dom2, num);
 #endif
 	  }
-	            
+
           MElement *e = createElementMSH(this, num, type, physical, elementary,
                                          partition, vertices, elements, physicals,
                                          own, p, doms[0], doms[1]);
-	  
+
 #if (FAST_ELEMENTS==1)
 	  elems[num] = e;
 	  elemreg[num] = elementary;
 	  elemphy[num] = physical;
-#endif	  
-	  
+#endif
+
           for(unsigned int j = 0; j < ghosts.size(); j++)
             _ghostCells.insert(std::pair<MElement*, short>(e, ghosts[j]));
           if(numElements > 100000)
             Msg::ProgressMeter(i + 1, numElements, "Reading elements");
           delete [] indices;
         }
-        
+
 #if (FAST_ELEMENTS==1)
 	for(int i = 0; i < 10; i++)
 	  elements[i].clear();
-	
+
 	std::map<int, MElement* >::iterator ite;
 	for (ite = elems.begin(); ite != elems.end(); ite++)
-	{	    
+	{
 	    int num = ite->first;
 	    MElement *e = ite->second;
 	    if (parents.find(num) == parents.end())
@@ -709,7 +709,7 @@ int GModel::readMSH(const std::string &name)
   // store the physical tags
   for(int i = 0; i < 4; i++)
     _storePhysicalTagsInEntities(i, physicals[i]);
-  
+
   fclose(fp);
 
   return postpro ? 2 : 1;
@@ -1514,7 +1514,7 @@ int GModel::readPLY(const std::string &name)
 
 	  pEnd = pEnd3;
           std::vector<double> prop(nbView);
-	  for (int k=0; k< nbView; k++){
+	  for (int k = 0; k < nbView; k++){
 	    prop[k]=strtod(pEnd, &pEnd2);
 	    pEnd = pEnd2;
 	    properties[k].push_back(prop[k]);
@@ -1524,7 +1524,7 @@ int GModel::readPLY(const std::string &name)
 	for(int i = 0; i < nbf; i++) {
 	  if(!fgets(buffer, sizeof(buffer), fp)) break;
 	  int n[3], nbe;
-	  int nb = sscanf(buffer, "%d %d %d %d", &nbe, &n[0], &n[1], &n[2]);
+	  sscanf(buffer, "%d %d %d %d", &nbe, &n[0], &n[1], &n[2]);
 	  std::vector<MVertex*> vertices;
 	  if(!getVertices(3, n, vertexVector, vertices)) return 0;
 	  elements[0][elementary].push_back(new MTriangle(vertices));
@@ -2324,7 +2324,7 @@ int GModel::writeMAIL(const std::string &name, bool saveAll, double scalingFacto
   for(fiter it = firstFace(); it != lastFace(); ++it){
     if(saveAll || (*it)->physicals.size()){
       for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
-        MTriangle *t = (*it)->triangles[i];
+        //MTriangle *t = (*it)->triangles[i];
         fprintf(fp, " %d %d %d\n", 0, 0, 0);
       }
     }
@@ -3190,7 +3190,6 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
       int v0, v1;
       char line[100000], *p, *pEnd, *pEnd2;
       int iLine = 1;
-      int num = 0;
       for (int k= 0; k < numElements; k++){
 	physicals[1][iLine][1] = "centerline";
 	if(!fgets(line, sizeof(line), fp)) return 0;
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index f746e6635f9bb8ceb2a35495d64d778b43c8c1d8..557113a6720a6332be4373b4c7fafa10347fc4a4 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -9,7 +9,7 @@
 #include "polynomialBasis.h"
 
 int MHexahedron::getVolumeSign()
-{ 
+{
   double mat[3][3];
   mat[0][0] = _v[1]->x() - _v[0]->x();
   mat[0][1] = _v[3]->x() - _v[0]->x();
@@ -40,7 +40,7 @@ void MHexahedron::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &r
     MVertex *v2 = _v[faces_hexa(ithFace, 2)];
     MVertex *v3 = _v[faces_hexa(ithFace, 3)];
 
-    if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && 
+    if (v0 == face.getVertex(0) && v1 == face.getVertex(1) &&
         v2 == face.getVertex(2) && v3 == face.getVertex(3)){
       sign = 1; rot = 0; return;
     }
@@ -109,7 +109,7 @@ void MHexahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3
   z[0] = pnt1.z(); z[1] = pnt2.z();
 
   // not great, but better than nothing
-  static const int f[6] = {0, 0, 0, 1, 2, 3};
+  //static const int f[6] = {0, 0, 0, 1, 2, 3};
   n[0] = n[1] = 1 ;
 }
 
@@ -196,116 +196,116 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl
   double oy = -1. + d*iy;
 
   if (io == 0){
-    double U1 = 
-      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25; 
-    double V1 = 
-      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25; 
-    double W1 = 
-      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25; 
+    double U1 =
+      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+    double V1 =
+      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+    double W1 =
+      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
 
     ox += d;
-  
-    double U2 = 
-      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25; 
-    double V2 = 
-      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25; 
-    double W2 = 
-      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25; 
+
+    double U2 =
+      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+    double V2 =
+      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+    double W2 =
+      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
 
     oy += d;
-  
-    double U3 = 
-      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25; 
-    double V3 = 
-      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25; 
-    double W3 = 
-      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25; 
-    
+
+    double U3 =
+      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+    double V3 =
+      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+    double W3 =
+      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
     hex->pnt(U1, V1, W1, pnt1);
     hex->pnt(U2, V2, W2, pnt2);
     hex->pnt(U3, V3, W3, pnt3);
-  } 
+  }
   else{
-    double U1 = 
-      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25; 
-    double V1 = 
-      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25; 
-    double W1 = 
-      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25; 
+    double U1 =
+      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+    double V1 =
+      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+    double W1 =
+      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
 
     ox += d;
     oy += d;
-  
-    double U2 = 
-      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25; 
-    double V2 = 
-      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25; 
-    double W2 = 
-      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25; 
+
+    double U2 =
+      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+    double V2 =
+      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+    double W2 =
+      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
 
     ox -= d;
-  
-    double U3 = 
-      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25; 
-    double V3 = 
-      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25; 
-    double W3 = 
-      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 + 
-      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 + 
-      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 + 
-      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25; 
-    
+
+    double U3 =
+      pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+    double V3 =
+      pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+    double W3 =
+      pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+      pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+      pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+      pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
     hex->pnt(U1, V1, W1, pnt1);
     hex->pnt(U2, V2, W2, pnt2);
     hex->pnt(U3, V3, W3, pnt3);
@@ -329,6 +329,6 @@ void MHexahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector3
 }
 
 int MHexahedronN::getNumFacesRep()
-{ 
+{
   return 6 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
 }