diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index 8ea4ce58f1bfb9014e6db2fff128d9362e37ccac..1345fbe1af71458481a9be1bf91add30d45be1d6 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -358,7 +358,7 @@ void CreateOutputFile(const std::string &fileName, int format, bool redraw)
     {
       if(!FlGui::available()) break;
 
-      FILE *fp = fopen(name.c_str(), "wb");
+      FILE *fp = Fopen(name.c_str(), "wb");
       if(!fp){
         Msg::Error("Unable to open file '%s'", name.c_str());
         error = true;
@@ -395,7 +395,7 @@ void CreateOutputFile(const std::string &fileName, int format, bool redraw)
     {
       if(!FlGui::available()) break;
 
-      FILE *fp = fopen(name.c_str(), "wb");
+      FILE *fp = Fopen(name.c_str(), "wb");
       if(!fp){
         Msg::Error("Unable to open file '%s'", name.c_str());
         error = true;
@@ -466,7 +466,7 @@ void CreateOutputFile(const std::string &fileName, int format, bool redraw)
     {
       if(!FlGui::available()) break;
 
-      FILE *fp = fopen(name.c_str(), "w");
+      FILE *fp = Fopen(name.c_str(), "w");
       if(!fp){
         Msg::Error("Unable to open file '%s'", name.c_str());
         error = true;
@@ -498,7 +498,7 @@ void CreateOutputFile(const std::string &fileName, int format, bool redraw)
   case FORMAT_MPEG:
     {
       std::string parFileName = CTX::instance()->homeDir + ".gmsh-mpeg_encode.par";
-      FILE *fp = fopen(parFileName.c_str(), "w");
+      FILE *fp = Fopen(parFileName.c_str(), "w");
       if(!fp){
         Msg::Error("Unable to open file '%s'", parFileName.c_str());
         error = true;
diff --git a/Common/GmshMessage.cpp b/Common/GmshMessage.cpp
index 10d2e18e61522029a2cc603927249fbffb72bb87..055a87cc73468d703054605fef6d54b2a865cf45 100644
--- a/Common/GmshMessage.cpp
+++ b/Common/GmshMessage.cpp
@@ -758,7 +758,7 @@ void Msg::InitializeOnelab(const std::string &name, const std::string &sockname)
   if(sockname.empty()){
     _onelabClient = new localGmsh();
     if(name != "Gmsh"){ // load db from file:
-      FILE *fp = fopen(name.c_str(), "rb");
+      FILE *fp = Fopen(name.c_str(), "rb");
       if(fp){
         _onelabClient->fromFile(fp);
         fclose(fp);
diff --git a/Common/OS.cpp b/Common/OS.cpp
index b7588ddf0711f0be55dab6ec3f290f980628bc08..9120f89fc99e8408e7fe3973fa5a8e9134e48b72 100644
--- a/Common/OS.cpp
+++ b/Common/OS.cpp
@@ -34,10 +34,6 @@
 #include <fstream>
 #endif
 
-#if defined(HAVE_FLTK)
-#include <FL/Fl.H> // for fl_fopen
-#endif
-
 #if defined(__APPLE__)
 #define RUSAGE_SELF      0
 #define RUSAGE_CHILDREN -1
@@ -45,11 +41,131 @@
 
 #include "GmshMessage.h"
 
-const char *GetEnvironmentVar(const char *var)
+#if defined(WIN32) && !defined(__CYGWIN__)
+
+// UTF8 utility routines borrowed from FLTK
+
+static unsigned utf8decode(const char* p, const char* end, int* len)
 {
-#if !defined(WIN32)
-  return getenv(var);
+  static unsigned short cp1252[32] = {
+    0x20ac, 0x0081, 0x201a, 0x0192, 0x201e, 0x2026, 0x2020, 0x2021,
+    0x02c6, 0x2030, 0x0160, 0x2039, 0x0152, 0x008d, 0x017d, 0x008f,
+    0x0090, 0x2018, 0x2019, 0x201c, 0x201d, 0x2022, 0x2013, 0x2014,
+    0x02dc, 0x2122, 0x0161, 0x203a, 0x0153, 0x009d, 0x017e, 0x0178
+  };
+  unsigned char c = *(unsigned char*)p;
+  if (c < 0x80) {
+    if (len) *len = 1;
+    return c;
+  } else if (c < 0xa0) {
+    if (len) *len = 1;
+    return cp1252[c-0x80];
+  } else if (c < 0xc2) {
+    goto FAIL;
+  }
+  if ( (end && p+1 >= end) || (p[1]&0xc0) != 0x80) goto FAIL;
+  if (c < 0xe0) {
+    if (len) *len = 2;
+    return
+      ((p[0] & 0x1f) << 6) +
+      ((p[1] & 0x3f));
+  } else if (c == 0xe0) {
+    if (((unsigned char*)p)[1] < 0xa0) goto FAIL;
+    goto UTF8_3;
+  } else if (c < 0xf0) {
+  UTF8_3:
+    if ( (end && p+2 >= end) || (p[2]&0xc0) != 0x80) goto FAIL;
+    if (len) *len = 3;
+    return
+      ((p[0] & 0x0f) << 12) +
+      ((p[1] & 0x3f) << 6) +
+      ((p[2] & 0x3f));
+  } else if (c == 0xf0) {
+    if (((unsigned char*)p)[1] < 0x90) goto FAIL;
+    goto UTF8_4;
+  } else if (c < 0xf4) {
+  UTF8_4:
+    if ( (end && p+3 >= end) || (p[2]&0xc0) != 0x80 || (p[3]&0xc0) != 0x80) goto FAIL;
+    if (len) *len = 4;
+    return
+      ((p[0] & 0x07) << 18) +
+      ((p[1] & 0x3f) << 12) +
+      ((p[2] & 0x3f) << 6) +
+      ((p[3] & 0x3f));
+  } else if (c == 0xf4) {
+    if (((unsigned char*)p)[1] > 0x8f) goto FAIL; /* after 0x10ffff */
+    goto UTF8_4;
+  } else {
+  FAIL:
+    if (len) *len = 1;
+    return c;
+  }
+}
+
+static unsigned utf8toUtf16(const char* src, unsigned srclen,
+                            unsigned short* dst, unsigned dstlen)
+{
+  const char* p = src;
+  const char* e = src+srclen;
+  unsigned count = 0;
+  if (dstlen) for (;;) {
+    if (p >= e) {dst[count] = 0; return count;}
+    if (!(*p & 0x80)) { /* ascii */
+      dst[count] = *p++;
+    } else {
+      int len; unsigned ucs = utf8decode(p,e,&len);
+      p += len;
+      if (ucs < 0x10000) {
+	dst[count] = ucs;
+      } else {
+	/* make a surrogate pair: */
+	if (count+2 >= dstlen) {dst[count] = 0; count += 2; break;}
+	dst[count] = (((ucs-0x10000u)>>10)&0x3ff) | 0xd800;
+	dst[++count] = (ucs&0x3ff) | 0xdc00;
+      }
+    }
+    if (++count == dstlen) {dst[count-1] = 0; break;}
+  }
+  /* we filled dst, measure the rest: */
+  while (p < e) {
+    if (!(*p & 0x80)) p++;
+    else {
+      int len; unsigned ucs = utf8decode(p,e,&len);
+      p += len;
+      if (ucs >= 0x10000) ++count;
+    }
+    ++count;
+  }
+  return count;
+}
+
+static wchar_t *wbuf = NULL;
+static wchar_t *wbuf1 = NULL;
+
+#endif
+
+FILE *Fopen(const char* f, const char *mode)
+{
+#if defined (WIN32) && !defined(__CYGWIN__)
+  size_t l = strlen(f);
+  unsigned wn = utf8toUtf16(f, (unsigned) l, NULL, 0) + 1;
+  wbuf = (wchar_t*)realloc(wbuf, sizeof(wchar_t)*wn);
+  wn = utf8toUtf16(f, (unsigned) l, (unsigned short *)wbuf, wn);
+  wbuf[wn] = 0;
+  l = strlen(mode);
+  wn = utf8toUtf16(mode, (unsigned) l, NULL, 0) + 1;
+  wbuf1 = (wchar_t*)realloc(wbuf1, sizeof(wchar_t)*wn);
+  wn = utf8toUtf16(mode, (unsigned) l, (unsigned short *)wbuf1, wn);
+  wbuf1[wn] = 0;
+  return _wfopen(wbuf, wbuf1);
 #else
+  return fopen(f, mode);
+#endif
+}
+
+const char *GetEnvironmentVar(const char *var)
+{
+#if defined(WIN32) && !defined(__CYGWIN__)
   const char *tmp = getenv(var);
   // Don't accept top dir or anything partially expanded like
   // c:\Documents and Settings\%USERPROFILE%, etc.
@@ -57,55 +173,48 @@ const char *GetEnvironmentVar(const char *var)
     return 0;
   else
     return tmp;
+#else
+  return getenv(var);
 #endif
 }
 
 const void SetEnvironmentVar(const char *var, const char *val)
 {
-#if !defined(WIN32)
-  setenv(var, val, 1);
-#else
+#if defined(WIN32) && !defined(__CYGWIN__)
   _putenv((std::string(var) + "=" + std::string(val)).c_str());
+#else
+  setenv(var, val, 1);
 #endif
 }
 
 double GetTimeInSeconds()
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
-  struct timeval tp;
-  gettimeofday(&tp, (struct timezone *)0);
-  double t = (double)tp.tv_sec + 1.e-6 * (double)tp.tv_usec;
-  return t;
-#else
+#if defined(WIN32) && !defined(__CYGWIN__)
   FILETIME ft;
   GetSystemTimeAsFileTime(&ft);
   double t =  1.e-7 * 4294967296. * (double)ft.dwHighDateTime +
               1.e-7 * (double)ft.dwLowDateTime;
   return t;
+#else
+  struct timeval tp;
+  gettimeofday(&tp, (struct timezone *)0);
+  double t = (double)tp.tv_sec + 1.e-6 * (double)tp.tv_usec;
+  return t;
 #endif
 }
 
 void SleepInSeconds(double s)
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
-  usleep((long)(1.e6 * s));
-#else
+#if defined(WIN32) && !defined(__CYGWIN__)
   Sleep((long)(1.e3 * s));
+#else
+  usleep((long)(1.e6 * s));
 #endif
 }
 
 static void GetResources(double *s, long *mem)
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
-  static struct rusage r;
-  getrusage(RUSAGE_SELF, &r);
-  *s = (double)r.ru_utime.tv_sec + 1.e-6 * (double)r.ru_utime.tv_usec;
-#if defined(__APPLE__)
-  *mem = (long)r.ru_maxrss;
-#else
-  *mem = (long)(r.ru_maxrss * 1024L);
-#endif
-#else
+#if defined(WIN32) && !defined(__CYGWIN__)
   FILETIME creation, exit, kernel, user;
   if(GetProcessTimes(GetCurrentProcess(), &creation, &exit, &kernel, &user)){
     *s = 1.e-7 * 4294967296. * (double)user.dwHighDateTime +
@@ -114,6 +223,15 @@ static void GetResources(double *s, long *mem)
   PROCESS_MEMORY_COUNTERS info;
   GetProcessMemoryInfo(GetCurrentProcess(), &info, sizeof(info));
   *mem = (long)info.PeakWorkingSetSize;
+#else
+  static struct rusage r;
+  getrusage(RUSAGE_SELF, &r);
+  *s = (double)r.ru_utime.tv_sec + 1.e-6 * (double)r.ru_utime.tv_usec;
+#if defined(__APPLE__)
+  *mem = (long)r.ru_maxrss;
+#else
+  *mem = (long)(r.ru_maxrss * 1024L);
+#endif
 #endif
 }
 
@@ -153,10 +271,10 @@ long GetMemoryUsage()
 
 int GetProcessId()
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
-  return getpid();
-#else
+#if defined(WIN32) && !defined(__CYGWIN__)
   return _getpid();
+#else
+  return getpid();
 #endif
 }
 
@@ -169,49 +287,61 @@ std::string GetHostName()
 
 int UnlinkFile(const std::string &fileName)
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
-  return unlink(fileName.c_str());
+#if defined(WIN32) && !defined(__CYGWIN__)
+  size_t l = strlen(fileName.c_str());
+  unsigned wn = utf8toUtf16(fileName.c_str(), (unsigned) l, NULL, 0) + 1;
+  wbuf = (wchar_t*)realloc(wbuf, sizeof(wchar_t)*wn);
+  wn = utf8toUtf16(fileName.c_str(), (unsigned) l, (unsigned short *)wbuf, wn);
+  wbuf[wn] = 0;
+  return _wunlink(wbuf);
 #else
-  return _unlink(fileName.c_str());
+  return unlink(fileName.c_str());
 #endif
 }
 
 int StatFile(const std::string &fileName)
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
+#if defined(WIN32) && !defined(__CYGWIN__)
+  struct _stat buf;
+  size_t l = strlen(fileName.c_str());
+  unsigned wn = utf8toUtf16(fileName.c_str(), (unsigned) l, NULL, 0) + 1;
+  wbuf = (wchar_t*)realloc(wbuf, sizeof(wchar_t)*wn);
+  wn = utf8toUtf16(fileName.c_str(), (unsigned) l, (unsigned short *)wbuf, wn);
+  wbuf[wn] = 0;
+  int ret = _wstat(wbuf, &buf);
+#else
   struct stat buf;
   int ret = stat(fileName.c_str(), &buf);
-  // could get file modification time from buf
-#else
-  struct _stat buf;
-  int ret = _stat(fileName.c_str(), &buf);
 #endif
   return ret;
 }
 
 int CreateDirectory(const std::string &dirName)
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
-  if(mkdir(dirName.c_str(), 0777))
-    return 0;
+#if defined(WIN32) && !defined(__CYGWIN__)
+  size_t l = strlen(dirName.c_str());
+  unsigned wn = utf8toUtf16(dirName.c_str(), (unsigned) l, NULL, 0) + 1;
+  wbuf = (wchar_t*)realloc(wbuf, sizeof(wchar_t)*wn);
+  wn = utf8toUtf16(dirName.c_str(), (unsigned) l, (unsigned short *)wbuf, wn);
+  wbuf[wn] = 0;
+  if(_wmkdir(wbuf)) return 0;
 #else
-  if(_mkdir(dirName.c_str()))
-    return 0;
+  if(mkdir(dirName.c_str(), 0777)) return 0;
 #endif
   return 1;
 }
 
 int KillProcess(int pid)
 {
-#if !defined(WIN32) || defined(__CYGWIN__)
-  if(kill(pid, 9))
-    return 0;
-#else
+#if defined(WIN32) && !defined(__CYGWIN__)
   HANDLE hProc = OpenProcess(PROCESS_TERMINATE, FALSE, pid);
   if(!TerminateProcess(hProc, 0)){
     CloseHandle(hProc);
     return 0;
   }
+#else
+  if(kill(pid, 9))
+    return 0;
 #endif
   return 1;
 }
@@ -242,7 +372,7 @@ int SystemCall(const std::string &command, bool blocking)
     }
   }
 
-#if defined(WIN32)
+#if defined(WIN32) && !defined(__CYGWIN__)
   if(isPython){
     Msg::Info("Shell opening '%s' with arguments '%s'", exe.c_str(),
               args.c_str());
@@ -303,11 +433,13 @@ int SystemCall(const std::string &command, bool blocking)
 std::string GetCurrentWorkdir()
 {
   char path[1024];
-#if defined(WIN32)
+
+#if defined(WIN32) && !defined(__CYGWIN__)
   if(!_getcwd(path, sizeof(path))) return "";
 #else
   if(!getcwd(path, sizeof(path))) return "";
 #endif
+
   std::string str(path);
   // match the convention of SplitFileName that delivers directory path
   // ending with a directory separator
@@ -372,13 +504,3 @@ void RedirectIOToConsole()
   std::ios::sync_with_stdio();
 #endif
 }
-
-FILE *Fopen(const char* f, const char *mode)
-{
-#if defined(HAVE_FLTK)
-  // this handles non-ASCII characters correctly on Windows
-  return fl_fopen(f, mode);
-#else
-  return fopen(f, mode);
-#endif
-}
diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index 9056ce296e3ee06b4d025caf29c4cfdf5d7e24ac..5e1ec88ffebfd295d9966578b56697400887988c 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -185,7 +185,7 @@ int ParseFile(const std::string &fileName, bool close, bool warnIfMissing)
   // add 'b' for pure Windows programs: opening in text mode messes up
   // fsetpos/fgetpos (used e.g. for user-defined functions)
   FILE *fp;
-  if(!(fp = fopen(fileName.c_str(), "rb"))){
+  if(!(fp = Fopen(fileName.c_str(), "rb"))){
     if(warnIfMissing)
       Msg::Warning("Unable to open file '%s'", fileName.c_str());
     return 0;
@@ -244,7 +244,7 @@ void ParseString(const std::string &str)
 {
   if(str.empty()) return;
   std::string fileName = CTX::instance()->homeDir + CTX::instance()->tmpFileName;
-  FILE *fp = fopen(fileName.c_str(), "w");
+  FILE *fp = Fopen(fileName.c_str(), "w");
   if(fp){
     fprintf(fp, "%s\n", str.c_str());
     fclose(fp);
@@ -276,7 +276,7 @@ int MergeFile(const std::string &fileName, bool warnIfMissing)
 
   // added 'b' for pure Windows programs, since some of these files
   // contain binary data
-  FILE *fp = fopen(fileName.c_str(), "rb");
+  FILE *fp = Fopen(fileName.c_str(), "rb");
   if(!fp){
     if(warnIfMissing)
       Msg::Warning("Unable to open file '%s'", fileName.c_str());
@@ -504,7 +504,7 @@ int MergePostProcessingFile(const std::string &fileName, bool showLastStep,
 {
 #if defined(HAVE_POST)
   // check if there is a mesh in the file
-  FILE *fp = fopen(fileName.c_str(), "rb");
+  FILE *fp = Fopen(fileName.c_str(), "rb");
   if(!fp){
     if(warnIfMissing) Msg::Warning("Unable to open file '%s'", fileName.c_str());
     return 0;
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 81e0c2556e83254369d20418fe3c3f5eacb73a23..1106ad88cfb1d9914e50f40d1ee5814ab9836bea 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -13,6 +13,7 @@
 #include "GModel.h"
 #include "Context.h"
 #include "Options.h"
+#include "OS.h"
 #include "Colors.h"
 #include "CommandLine.h"
 #include "GamePad.h"
@@ -617,7 +618,7 @@ void PrintOptions(int num, int level, int diff, int help, const char *filename,
   FILE *file;
 
   if(filename) {
-    file = fopen(filename, "w");
+    file = Fopen(filename, "w");
     if(!file) {
       Msg::Error("Unable to open file '%s'", filename);
       return;
@@ -732,7 +733,7 @@ void PrintOptionsDoc()
     "@c\n\n";
 
   {
-    FILE *file = fopen("opt_general.texi", "w");
+    FILE *file = Fopen("opt_general.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_general.texi'");
       return;
@@ -745,7 +746,7 @@ void PrintOptionsDoc()
     fclose(file);
   }
   {
-    FILE *file = fopen("opt_print.texi", "w");
+    FILE *file = Fopen("opt_print.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_print.texi'");
       return;
@@ -758,7 +759,7 @@ void PrintOptionsDoc()
     fclose(file);
   }
   {
-    FILE *file = fopen("opt_geometry.texi", "w");
+    FILE *file = Fopen("opt_geometry.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_geometry.texi'");
       return;
@@ -771,7 +772,7 @@ void PrintOptionsDoc()
     fclose(file);
   }
   {
-    FILE *file = fopen("opt_mesh.texi", "w");
+    FILE *file = Fopen("opt_mesh.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_mesh.texi'");
       return;
@@ -784,7 +785,7 @@ void PrintOptionsDoc()
     fclose(file);
   }
   {
-    FILE *file = fopen("opt_solver.texi", "w");
+    FILE *file = Fopen("opt_solver.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_solver.texi'");
       return;
@@ -797,7 +798,7 @@ void PrintOptionsDoc()
     fclose(file);
   }
   {
-    FILE *file = fopen("opt_post.texi", "w");
+    FILE *file = Fopen("opt_post.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_post.texi'");
       return;
@@ -811,7 +812,7 @@ void PrintOptionsDoc()
   }
   {
 #if defined(HAVE_POST)
-    FILE *file = fopen("opt_view.texi", "w");
+    FILE *file = Fopen("opt_view.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_view.texi'");
       return;
@@ -830,7 +831,7 @@ void PrintOptionsDoc()
   }
   {
 #if defined(HAVE_PLUGINS)
-    FILE *file = fopen("opt_plugin.texi", "w");
+    FILE *file = Fopen("opt_plugin.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_plugin.texi'");
       return;
@@ -874,7 +875,7 @@ void PrintOptionsDoc()
 
 #if defined(HAVE_MESH)
   {
-    FILE *file = fopen("opt_fields.texi", "w");
+    FILE *file = Fopen("opt_fields.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'opt_fields.texi'");
       return;
@@ -919,7 +920,7 @@ void PrintOptionsDoc()
   }
 #endif
   {
-    FILE *file = fopen("shortcuts.texi", "w");
+    FILE *file = Fopen("shortcuts.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'shortcuts.texi'");
       return;
@@ -932,7 +933,7 @@ void PrintOptionsDoc()
     fclose(file);
   }
   {
-    FILE *file = fopen("mouse.texi", "w");
+    FILE *file = Fopen("mouse.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'mouse.texi'");
       return;
@@ -945,7 +946,7 @@ void PrintOptionsDoc()
     fclose(file);
   }
   {
-    FILE *file = fopen("commandline.texi", "w");
+    FILE *file = Fopen("commandline.texi", "w");
     if(!file) {
       Msg::Error("Unable to open file 'commandline.texi'");
       return;
diff --git a/Common/SmoothData.cpp b/Common/SmoothData.cpp
index 3ca7e75cb7c6edeef05739614be087cf8cabbb09..c81bbc5d7bd13611ff1bb00fe800db255896a7cc 100644
--- a/Common/SmoothData.cpp
+++ b/Common/SmoothData.cpp
@@ -6,13 +6,14 @@
 #include <stdio.h>
 #include "SmoothData.h"
 #include "Numeric.h"
+#include "OS.h"
 
 // Basic coordinate-based floting point value averager
 
 double xyzv::eps = 1.e-12;
 
 xyzv::xyzv(const xyzv &other)
-{ 
+{
   x = other.x;
   y = other.y;
   z = other.z;
@@ -42,7 +43,7 @@ xyzv &xyzv::operator = (const xyzv &other)
   return *this;
 }
 
-void xyzv::update(int n, double *v) 
+void xyzv::update(int n, double *v)
 {
   if(!vals) {
     vals = new double[n];
@@ -97,14 +98,14 @@ void smooth_data::normalize()
 
 bool smooth_data::exportview(std::string filename)
 {
-  FILE *fp = fopen(filename.c_str(), "w");
+  FILE *fp = Fopen(filename.c_str(), "w");
   if(!fp) return false;
   fprintf(fp, "View \"data\" {\n");
   std::set<xyzv, lessthanxyzv>::iterator it = c.begin();
   while(it != c.end()){
     switch(it->nbvals){
-    case 1: 
-      fprintf(fp, "SP(%.16g,%.16g,%.16g){%.16g};\n", 
+    case 1:
+      fprintf(fp, "SP(%.16g,%.16g,%.16g){%.16g};\n",
               it->x, it->y, it->z, it->vals[0]);
       break;
     case 3:
@@ -127,17 +128,17 @@ float xyzn::angle(int i, char nx, char ny, char nz)
 {
   // returns the angle (in [-180,180]) between the ith normal stored
   // at point xyz and the new normal nx,ny,nz
-  double a[3] = {char2float(n[i].nx), 
-                 char2float(n[i].ny), 
+  double a[3] = {char2float(n[i].nx),
+                 char2float(n[i].ny),
                  char2float(n[i].nz)};
-  double b[3] = {char2float(nx), 
-                 char2float(ny), 
+  double b[3] = {char2float(nx),
+                 char2float(ny),
                  char2float(nz)};
   norme(a);
   norme(b);
   double c[3];
   prodve(a, b, c);
-  double cosc; 
+  double cosc;
   prosca(a, b, &cosc);
   double sinc = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
   double angplan = myatan2(sinc, cosc);
@@ -181,7 +182,7 @@ void smooth_normals::add(double x, double y, double z,
   xyzn xyz((float)x, (float)y, (float)z);
   std::set<xyzn, lessthanxyzn>::const_iterator it = c.find(xyz);
   if(it == c.end()) {
-    xyz.update(float2char((float)nx), 
+    xyz.update(float2char((float)nx),
                float2char((float)ny),
                float2char((float)nz), tol);
     c.insert(xyz);
@@ -190,23 +191,23 @@ void smooth_normals::add(double x, double y, double z,
     // we can do this because we know that it will not destroy the set
     // ordering
     xyzn *p = (xyzn *) & (*it);
-    p->update(float2char((float)nx), 
-              float2char((float)ny), 
+    p->update(float2char((float)nx),
+              float2char((float)ny),
               float2char((float)nz), tol);
-  }    
+  }
 }
 
 bool smooth_normals::get(double x, double y, double z,
                          double &nx, double &ny, double &nz)
 {
-  std::set<xyzn, lessthanxyzn>::const_iterator it = 
+  std::set<xyzn, lessthanxyzn>::const_iterator it =
     c.find(xyzn((float)x, (float)y, (float)z));
 
   if(it == c.end()) return false;
 
   xyzn *p = (xyzn *) & (*it);
   for(unsigned int i = 0; i < p->n.size(); i++){
-    if(fabs(p->angle(i, float2char((float)nx), 
+    if(fabs(p->angle(i, float2char((float)nx),
                      float2char((float)ny),
                      float2char((float)nz))) < tol) {
       nx = char2float(p->n[i].nx);
diff --git a/Common/StringUtils.cpp b/Common/StringUtils.cpp
index 22ada724fe850c67a01bc48c08fb808806375ae5..1b496dfbc77b77fd1f3f85246349acde79ee8a9a 100644
--- a/Common/StringUtils.cpp
+++ b/Common/StringUtils.cpp
@@ -7,9 +7,9 @@
 #if defined(__CYGWIN__)
 #include <sys/cygwin.h>
 #endif
-
 #include "StringUtils.h"
 #include "GmshMessage.h"
+#include "OS.h"
 
 void SwapBytes(char *array, int size, int n)
 {
@@ -112,7 +112,7 @@ std::string GetFileNameWithoutPath(const std::string &fileName)
 
 std::string ConvertFileToString(const std::string &fileName)
 {
-  FILE *fp = fopen(fileName.c_str(), "r");
+  FILE *fp = Fopen(fileName.c_str(), "r");
   if(!fp) return "";
   std::string out;
   char str[256];
@@ -121,7 +121,7 @@ std::string ConvertFileToString(const std::string &fileName)
   return out;
 }
 
-void ReplaceSubStringInPlace(const std::string &olds, const std::string &news, 
+void ReplaceSubStringInPlace(const std::string &olds, const std::string &news,
                              std::string &str)
 {
   while(1){
@@ -131,7 +131,7 @@ void ReplaceSubStringInPlace(const std::string &olds, const std::string &news,
   }
 }
 
-std::string ReplaceSubString(const std::string &olds, const std::string &news, 
+std::string ReplaceSubString(const std::string &olds, const std::string &news,
                              const std::string &str)
 {
   std::string copy(str);
diff --git a/Common/onelab.h b/Common/onelab.h
index f093c8bbf0b9ea9fab46c6ec52fa6a4c08814893..b99af7787429560398eee6d0db6ccb6269f45a9f 100644
--- a/Common/onelab.h
+++ b/Common/onelab.h
@@ -929,6 +929,7 @@ namespace onelab{
       if(!_server) _server = new server(address);
       return _server;
     }
+    static void setInstance(server *s) { _server = s; }
     void clear(const std::string &name="", const std::string &client="")
     {
       _parameterSpace.clear(name, client);
diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index a9bd8a29436d77014f1d305d138fc2ded751c332..d6ada68fc30df33ef9891b5f92345d577e16f36a 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -71,7 +71,7 @@ static void file_new_cb(Fl_Widget *w, void *data)
       else
         goto test;
     }
-    FILE *fp = fopen(name.c_str(), "w");
+    FILE *fp = Fopen(name.c_str(), "w");
     if(!fp){
       Msg::Error("Unable to open file '%s'", name.c_str());
       return;
@@ -3061,7 +3061,7 @@ void graphicWindow::saveMessages(const char *filename)
 {
   if(!_browser) return;
 
-  FILE *fp = fopen(filename, "w");
+  FILE *fp = Fopen(filename, "w");
 
   if(!fp) {
     Msg::Error("Unable to open file '%s'", filename);
diff --git a/Fltk/onelabGroup.cpp b/Fltk/onelabGroup.cpp
index 53ae44c7a936f425f5eff84f86041b2ff3907390..da18acb116d3521da2beaae26cc3d7c5e1820169 100644
--- a/Fltk/onelabGroup.cpp
+++ b/Fltk/onelabGroup.cpp
@@ -589,7 +589,7 @@ static std::string timeStamp()
 
 static void saveDb(const std::string &fileName)
 {
-  FILE *fp = fopen(fileName.c_str(), "wb");
+  FILE *fp = Fopen(fileName.c_str(), "wb");
   if(fp){
     Msg::StatusBar(true, "Saving database '%s'...", fileName.c_str());
     onelab::server::instance()->toFile(fp);
@@ -644,7 +644,7 @@ static void archiveOutputFiles(const std::string &fileName)
 static void loadDb(const std::string &name)
 {
   Msg::StatusBar(true, "Loading database '%s'...", name.c_str());
-  FILE *fp = fopen(name.c_str(), "rb");
+  FILE *fp = Fopen(name.c_str(), "rb");
   if(fp){
     onelab::server::instance()->fromFile(fp);
     fclose(fp);
diff --git a/Fltk/pluginWindow.cpp b/Fltk/pluginWindow.cpp
index 63a213e03f9334dcb9c5765453b46e8d38601a69..e36291f1c9fb5921f5542e4a4fda9cf55e554a04 100644
--- a/Fltk/pluginWindow.cpp
+++ b/Fltk/pluginWindow.cpp
@@ -31,6 +31,7 @@ typedef unsigned long intptr_t;
 #include "Context.h"
 #include "GeoStringInterface.h"
 #include "StringUtils.h"
+#include "OS.h"
 
 #define MAX_PLUGIN_OPTIONS 50
 class PluginDialogBox{
@@ -132,7 +133,7 @@ static void add_scripting(GMSH_PostPlugin *p, PView *view)
     fileName = GModel::current()->getFileName();
 
   fileName += ".opt";
-  FILE *fp = fopen(fileName.c_str(), "a");
+  FILE *fp = Fopen(fileName.c_str(), "a");
   if(!fp){
     Msg::Error("Could not open file '%s'", fileName.c_str());
   }
diff --git a/Fltk/projectionEditor.cpp b/Fltk/projectionEditor.cpp
index bf09f357786927aa5ee60a9e18205dd1bc8e2812..bac9260ee4af6208c268d83fa40ae0fc9ba51f7e 100644
--- a/Fltk/projectionEditor.cpp
+++ b/Fltk/projectionEditor.cpp
@@ -21,6 +21,7 @@
 #include "fourierFace.h"
 #include "GmshMessage.h"
 #include "Context.h"
+#include "OS.h"
 
 #if defined(HAVE_FOURIER_MODEL)
 
@@ -419,7 +420,7 @@ static void save_selection_cb(Fl_Widget *w, void *data)
 {
   projectionEditor *e = (projectionEditor*)data;
   if(fileChooser(FILE_CHOOSER_CREATE, "Save Selection", "*.{geo,msh}")){
-    FILE *fp = fopen(fileChooserGetName(1).c_str(), "w");
+    FILE *fp = Fopen(fileChooserGetName(1).c_str(), "w");
     if(!fp){
       Msg::Error("Unable to open file `%s'", fileChooserGetName(1).c_str());
       return;
@@ -459,7 +460,7 @@ static void load_projection_cb(Fl_Widget *w, void *data)
 {
   projectionEditor *e = (projectionEditor*)data;
   if(fileChooser(FILE_CHOOSER_SINGLE, "Load Projection", "*.pro")){
-    FILE *fp = fopen(fileChooserGetName(1).c_str(), "r");
+    FILE *fp = Fopen(fileChooserGetName(1).c_str(), "r");
     if(!fp){
       Msg::Error("Unable to open file `%s'", fileChooserGetName(1).c_str());
       return;
@@ -506,7 +507,7 @@ static void save_projection_cb(Fl_Widget *w, void *data)
     FM::ProjectionSurface *ps = (FM::ProjectionSurface*)p->face->getNativePtr();
     if(fileChooser(FILE_CHOOSER_CREATE, "Save Projection", "*.pro")){
       std::string name = fileChooserGetName(1);
-      FILE *fp = fopen(name.c_str(), "w");
+      FILE *fp = Fopen(name.c_str(), "w");
       if(!fp){
         Msg::Error("Unable to open file `%s'", name.c_str());
         return;
@@ -691,7 +692,7 @@ static void action_cb(Fl_Widget *w, void *data)
   }
   else{
     if(fileChooser(FILE_CHOOSER_CREATE, "Save Fourier Model", "*.fm")){
-      FILE *fp = fopen(fileChooserGetName(1).c_str(), "w");
+      FILE *fp = Fopen(fileChooserGetName(1).c_str(), "w");
       if(!fp){
         Msg::Error("Unable to open file `%s'", fileChooserGetName(1).c_str());
         return;
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index a2dda55db7db6aeef6006c5a89d1bfe9f3128ed9..4b978a81c32fbd64402f907b453dd77a2a0d08d6 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -269,7 +269,7 @@ void statisticsWindow::compute(bool elementQuality)
   //   int nbEdges = edges.size();
   //   printf("nb edges =%d \n", nbEdges);
   //   if(system("rm qualEdges.txt"));
-  //   FILE *fp = fopen("qualEdges.txt", "w");
+  //   FILE *fp = Fopen("qualEdges.txt", "w");
   //   std::vector<int> qualE;
   //   int nbS = 50;
   //   qualE.resize(nbS);
@@ -347,7 +347,7 @@ void statisticsWindow::compute(bool elementQuality)
   //   Field *f = fields->get(fields.getBackgroundField());
   //   int nbEdges = edges.size();
   //   if(system("rm qualEdges.txt"));
-  //   FILE *fp = fopen("qualEdges.txt", "w");
+  //   FILE *fp = Fopen("qualEdges.txt", "w");
   //   std::vector<int> qualE;
   //   int nbS = 50;
   //   qualE.resize(nbS);
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 43fe44cdd134200fa898cb1d32156800187001dc..f003cb68915659f06d3bbe1fadf8dbe9c015adea 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -967,7 +967,7 @@ void CellComplex::printComplex(int dim)
 
 int CellComplex::saveComplex(std::string filename)
 {
-  /*FILE *fp = fopen (filename.c_str(), "w");
+  /*FILE *fp = Fopen (filename.c_str(), "w");
   if(!fp){
     printf("\nUnable to open file '%s' \n", filename.c_str());
     return 0;
@@ -1008,7 +1008,7 @@ int CellComplex::saveComplex(std::string filename)
 
 int CellComplex::loadComplex(std::string filename)
 {
-  /*  FILE *fp = fopen (filename.c_str(), "r");
+  /*  FILE *fp = Fopen (filename.c_str(), "r");
   if(!fp){
     printf("\nUnable to open file '%s' \n", filename.c_str());
     return 0;
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index b848eb3eed548acfdd7ce15d73f3e13a188915b5..c64f8b000e8a8714ea3c4d79e96e369b7edceda2 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -17,6 +17,7 @@
 #include "Numeric.h"
 #include "GaussLegendre1D.h"
 #include "Context.h"
+#include "OS.h"
 
 #if defined(HAVE_MESH)
 #include "meshGFaceOptimize.h"
@@ -966,7 +967,7 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
     printf("dist = %12.5E\n",dist);
   }
 
-  // FILE *F = fopen ("hop.pos","w");
+  // FILE *F = Fopen ("hop.pos","w");
   // fprintf(F,"View \" \" {\n");
   // fprintf(F,"SP(%g,%g,%g) {%g};\n",queryPoint.x(),queryPoint.y(),queryPoint.z(),0.0);
   double initial_guesses = 10.0;
@@ -1406,7 +1407,7 @@ void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio
 {
   SVector3 ez (0, 0, 1);
   std::list<GEdgeLoop>::iterator it = edgeLoops.begin();
-  FILE *f = fopen ("coucou.pos","w");
+  FILE *f = Fopen ("coucou.pos","w");
   fprintf(f,"View \"\"{\n");
   for (; it != edgeLoops.end(); ++it){
     bool found = false;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index e68ae17b5959a96fc4c269baffbd7e72581d6480..078e119e9d994816f34521816862c7856ce2a6e9 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -12,6 +12,7 @@
 #include "GFaceCompound.h"
 #include "GEdgeCompound.h"
 #include "intersectCurveSurface.h"
+#include "OS.h"
 
 #if defined(HAVE_SOLVER) && defined(HAVE_ANN)
 
@@ -19,7 +20,6 @@
 #include "MLine.h"
 #include "MTriangle.h"
 #include "Numeric.h"
-#include "OS.h"
 #include "Octree.h"
 #include "SBoundingBox3d.h"
 #include "SPoint3.h"
@@ -61,7 +61,7 @@ static void printBound(std::vector<MVertex*> &l, int tag)
 {
   char name[256];
   sprintf(name, "myBOUND%d.pos", tag);
-  FILE * xyz = fopen(name,"w");
+  FILE * xyz = Fopen(name,"w");
   fprintf(xyz,"View \"\"{\n");
   for(unsigned int i = 0; i < l.size(); i++){
     MVertex *v = l[i];
@@ -454,7 +454,7 @@ void GFaceCompound::printFillTris() const
     char name[256];
     //std::list<GFace*>::const_iterator itf = _compound.begin();
     sprintf(name, "fillTris-%d.pos", tag());
-    FILE * ftri = fopen(name,"w");
+    FILE * ftri = Fopen(name,"w");
     fprintf(ftri,"View \"\"{\n");
     for (std::list<MTriangle*>::iterator it2 = fillTris.begin();
 	 it2 !=fillTris.end(); it2++ ){
@@ -1289,7 +1289,7 @@ GFaceCompound::~GFaceCompound()
   if(kdtree){
     ANNpointArray nodes = kdtree->thePoints();
     if(nodes) annDeallocPts(nodes);
-    delete kdtree; 
+    delete kdtree;
   }
 }
 
@@ -2802,12 +2802,12 @@ void GFaceCompound::printStuff(int iNewton) const
 
   sprintf(name7, "UVM-%d.pos", (*it)->tag());
 
-  //FILE * uva = fopen(name0,"w");
-  FILE * uvx = fopen(name1,"w");
-  FILE * uvy = fopen(name2,"w");
-  FILE * uvz = fopen(name3,"w");
-  FILE * xyzu = fopen(name4,"w");
-  FILE * xyzv = fopen(name5,"w");
+  //FILE * uva = Fopen(name0,"w");
+  FILE * uvx = Fopen(name1,"w");
+  FILE * uvy = Fopen(name2,"w");
+  FILE * uvz = Fopen(name3,"w");
+  FILE * xyzu = Fopen(name4,"w");
+  FILE * xyzv = Fopen(name5,"w");
 
   fprintf(uvx, "View \"\"{\n");
   fprintf(uvy, "View \"\"{\n");
diff --git a/Geo/GModelIO_ACIS.cpp b/Geo/GModelIO_ACIS.cpp
index d5c339ec88c0411c4d1c91bceebaed335adcfa31..fa3fc4bdbc6cf0a1768786097855f912b81432d0 100644
--- a/Geo/GModelIO_ACIS.cpp
+++ b/Geo/GModelIO_ACIS.cpp
@@ -6,6 +6,7 @@
 #include "GmshConfig.h"
 #include "GModel.h"
 #include "GmshMessage.h"
+#include "OS.h"
 #include "GModelIO_ACIS.h"
 #include "ACISVertex.h"
 #include "ACISEdge.h"
@@ -207,7 +208,7 @@ void ACIS_Internals::addFaces (GModel *gm, ENTITY_LIST &l)
 
 void ACIS_Internals::loadSAT(std::string fileName, GModel *gm)
 {
-  FILE *f = fopen (fileName.c_str(), "r");
+  FILE *f = Fopen (fileName.c_str(), "r");
   if (!f){
     return;
   }
diff --git a/Geo/GModelIO_BDF.cpp b/Geo/GModelIO_BDF.cpp
index e1c0a54d1a4107dcef42c073a097d22a12a76cbc..a2f07f40c7b18c6226476dcc776b27a927593278 100644
--- a/Geo/GModelIO_BDF.cpp
+++ b/Geo/GModelIO_BDF.cpp
@@ -6,6 +6,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include "GModel.h"
+#include "OS.h"
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
@@ -174,7 +175,7 @@ static int readElementBDF(FILE *fp, char *buffer, int keySize, int numVertices,
 
 int GModel::readBDF(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -301,7 +302,7 @@ int GModel::readBDF(const std::string &name)
 int GModel::writeBDF(const std::string &name, int format, int elementTagType,
                      bool saveAll, double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_DIFF.cpp b/Geo/GModelIO_DIFF.cpp
index 2ad72e9d0f24bfea52db9addc65a553b35f307d0..83b5dc511555b7c3cb0c110dcfd0952b45f82bfb 100644
--- a/Geo/GModelIO_DIFF.cpp
+++ b/Geo/GModelIO_DIFF.cpp
@@ -6,6 +6,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include "GModel.h"
+#include "OS.h"
 #include "MElement.h"
 
 static bool getVertices(int num, int *indices, std::map<int, MVertex*> &map,
@@ -38,7 +39,7 @@ static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
 
 int GModel::readDIFF(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -366,7 +367,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
     return 0;
   }
 
-  FILE *fp = fopen(name.c_str(), binary ? "wb" : "w");
+  FILE *fp = Fopen(name.c_str(), binary ? "wb" : "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index 97c6c9746dcf540027992981621b1bfe8d97c92f..2293dd0073fbd75af8e35982f5a5f45f98556cdf 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -8,6 +8,7 @@
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "GModel.h"
+#include "OS.h"
 #include "Geo.h"
 #include "OpenFile.h"
 #include "Numeric.h"
@@ -162,7 +163,7 @@ int GModel::importGEOInternals()
           e = new gmshEdge(this, c, 0, 0);
           add(e);
         }
-  
+
 
         if(!c->Visible) e->setVisibility(0);
         if(c->Color.type) e->setColor(c->Color.mesh);
@@ -452,7 +453,7 @@ static bool skipVertex(GVertex *gv)
 
 int GModel::writeGEO(const std::string &name, bool printLabels, bool onlyPhysicals)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
 
   std::map<double, std::string> meshSizeParameters;
   int cpt = 0;
diff --git a/Geo/GModelIO_GEOM.cpp b/Geo/GModelIO_GEOM.cpp
index 813ddc44aeb0827b2fb8da49efa076b5df3d3c33..e318fc8dd0d75d13a2b54c70051f0c0f124cd5b9 100644
--- a/Geo/GModelIO_GEOM.cpp
+++ b/Geo/GModelIO_GEOM.cpp
@@ -4,6 +4,7 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "OS.h"
 #include "MTriangle.h"
 
 static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
@@ -24,7 +25,7 @@ int GModel::readGEOM(const std::string &name)
 {
   // this is a format (from geomview?) that Bruno Levy's Graphite code
   // can write
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_INP.cpp b/Geo/GModelIO_INP.cpp
index 5f525aeb1de599151a7e7902b5e9d1a8325054b2..0c82a562aa687e30af318ec065b4bd3d2d5e962e 100644
--- a/Geo/GModelIO_INP.cpp
+++ b/Geo/GModelIO_INP.cpp
@@ -4,6 +4,7 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "OS.h"
 #include "MPoint.h"
 #include "MLine.h"
 #include "MTriangle.h"
@@ -46,7 +47,7 @@ static std::string physicalName(GModel *m, int dim, int num)
 int GModel::writeINP(const std::string &name, bool saveAll, bool saveGroupsOfNodes,
                      double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_IR3.cpp b/Geo/GModelIO_IR3.cpp
index 4b5b3ebb827cb76bf82e0a39ad2d950c59b5c269..2cc29ce9c58fbc3ce3fd581313a68d39c9f7afd2 100644
--- a/Geo/GModelIO_IR3.cpp
+++ b/Geo/GModelIO_IR3.cpp
@@ -4,12 +4,13 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "OS.h"
 #include "MElement.h"
 
 int GModel::writeIR3(const std::string &name, int elementTagType,
                      bool saveAll, double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_MAIL.cpp b/Geo/GModelIO_MAIL.cpp
index b6a2041639b4b95aa09b503107fd00342ac57134..1dda663233a7e9db4e4f80b953aa3eb5126bd496 100644
--- a/Geo/GModelIO_MAIL.cpp
+++ b/Geo/GModelIO_MAIL.cpp
@@ -4,6 +4,7 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "OS.h"
 #include "MTriangle.h"
 
 int GModel::writeMAIL(const std::string &name, bool saveAll, double scalingFactor)
@@ -11,7 +12,7 @@ int GModel::writeMAIL(const std::string &name, bool saveAll, double scalingFacto
   // CEA triangulation (.mail format) for Eric Darrigrand. Note that
   // we currently don't save the edges of the triangulation (the last
   // part of the file).
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_MESH.cpp b/Geo/GModelIO_MESH.cpp
index 49bc0fba00be526078bb21e5f868f19fcbf826a8..f23fd8b8970f75c8bd8e34f43ea4ba63dd62e250 100644
--- a/Geo/GModelIO_MESH.cpp
+++ b/Geo/GModelIO_MESH.cpp
@@ -6,6 +6,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include "GModel.h"
+#include "OS.h"
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
@@ -29,7 +30,7 @@ static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
 
 int GModel::readMESH(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -210,7 +211,7 @@ int GModel::readMESH(const std::string &name)
 int GModel::writeMESH(const std::string &name, int elementTagType,
                       bool saveAll, double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index 32fdcf9cedad25f987c844ab513fa6580cbbfd8b..764661a08b297c70c6b9585262f647a85373b22e 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/Geo/GModelIO_MSH.cpp
@@ -4,6 +4,7 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "OS.h"
 #include "GmshMessage.h"
 #include "MElement.h"
 #include "MPoint.h"
@@ -78,7 +79,7 @@ void readMSHPeriodicNodes(FILE *fp, GModel *gm)
 
 int GModel::readMSH(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "rb");
+  FILE *fp = Fopen(name.c_str(), "rb");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -472,9 +473,9 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
 
   FILE *fp;
   if(multipleView)
-    fp = fopen(name.c_str(), binary ? "ab" : "a");
+    fp = Fopen(name.c_str(), binary ? "ab" : "a");
   else
-    fp = fopen(name.c_str(), binary ? "wb" : "w");
+    fp = Fopen(name.c_str(), binary ? "wb" : "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_MSH2.cpp b/Geo/GModelIO_MSH2.cpp
index d3e1485d4bf2114c865e35f9100b21c26214df99..a082dc114f6a4d55c42fa6f96e4e64845b280e58 100644
--- a/Geo/GModelIO_MSH2.cpp
+++ b/Geo/GModelIO_MSH2.cpp
@@ -11,6 +11,7 @@
 #include <cassert>
 #include <iomanip>
 #include "GModel.h"
+#include "OS.h"
 #include "GmshDefines.h"
 #include "MPoint.h"
 #include "MLine.h"
@@ -180,7 +181,7 @@ static void getDomains(int dom1Num, int dom2Num, int type,
 
 int GModel::_readMSH2(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "rb");
+  FILE *fp = Fopen(name.c_str(), "rb");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -842,9 +843,9 @@ int GModel::_writeMSH2(const std::string &name, double version, bool binary,
 {
   FILE *fp;
   if(multipleView)
-    fp = fopen(name.c_str(), binary ? "ab" : "a");
+    fp = Fopen(name.c_str(), binary ? "ab" : "a");
   else
-    fp = fopen(name.c_str(), binary ? "wb" : "w");
+    fp = Fopen(name.c_str(), binary ? "wb" : "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -1069,7 +1070,7 @@ int GModel::writePartitionedMSH(const std::string &baseName, bool binary,
 #if 0
   if(_ghostCells.size()){
     Msg::Info("Writing ghost cells in debug file 'ghosts.pos'");
-    FILE *fp = fopen("ghosts.pos", "w");
+    FILE *fp = Fopen("ghosts.pos", "w");
     fprintf(fp, "View \"ghosts\"{\n");
     for(std::multimap<MElement*, short>::iterator it = _ghostCells.begin();
         it != _ghostCells.end(); it++)
diff --git a/Geo/GModelIO_P3D.cpp b/Geo/GModelIO_P3D.cpp
index eaa2947ea3da2012675b783c8071f767354d55dd..38ce548a3e7f5cff5951fe10202fdf929df002bf 100644
--- a/Geo/GModelIO_P3D.cpp
+++ b/Geo/GModelIO_P3D.cpp
@@ -4,6 +4,7 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "OS.h"
 #include "MQuadrangle.h"
 #include "MHexahedron.h"
 #include "discreteFace.h"
@@ -11,7 +12,7 @@
 
 int GModel::readP3D(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -110,7 +111,7 @@ int GModel::readP3D(const std::string &name)
 
 int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_PLY.cpp b/Geo/GModelIO_PLY.cpp
index b172b0fd5aa08d0bc1fa1268bb10c504f00372d8..949b419558dccaf664990f493416e9f9df15a898 100644
--- a/Geo/GModelIO_PLY.cpp
+++ b/Geo/GModelIO_PLY.cpp
@@ -57,7 +57,7 @@ int GModel::readPLY(const std::string &name)
   // this is crazy!?
   replaceCommaByDot(name);
 
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -177,7 +177,7 @@ int GModel::readPLY(const std::string &name)
 
 int GModel::readPLY2(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -240,7 +240,7 @@ int GModel::readPLY2(const std::string &name)
 
 int GModel::writePLY2(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_POS.cpp b/Geo/GModelIO_POS.cpp
index e58e1ff2fac5d2498b323911e3304b938aea59f9..0b0d3cbc4f0b50d5598bb3dfe7c275c1f23a6f8f 100644
--- a/Geo/GModelIO_POS.cpp
+++ b/Geo/GModelIO_POS.cpp
@@ -5,6 +5,7 @@
 
 #include <stdio.h>
 #include "GModel.h"
+#include "OS.h"
 #include "MElement.h"
 
 int GModel::writePOS(const std::string &name, bool printElementary,
@@ -12,7 +13,7 @@ int GModel::writePOS(const std::string &name, bool printElementary,
                      bool printRho, bool printDisto,
                      bool saveAll, double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_STL.cpp b/Geo/GModelIO_STL.cpp
index 4a7fb22430ce736362bfe954036de859d29ca6c7..9b0390cd3a3484e42dade937b5b505b5cd0a79a6 100644
--- a/Geo/GModelIO_STL.cpp
+++ b/Geo/GModelIO_STL.cpp
@@ -5,6 +5,7 @@
 
 #include <stdio.h>
 #include "GModel.h"
+#include "OS.h"
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
@@ -14,7 +15,7 @@
 
 int GModel::readSTL(const std::string &name, double tolerance)
 {
-  FILE *fp = fopen(name.c_str(), "rb");
+  FILE *fp = Fopen(name.c_str(), "rb");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -176,7 +177,7 @@ int GModel::readSTL(const std::string &name, double tolerance)
 int GModel::writeSTL(const std::string &name, bool binary, bool saveAll,
                      double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), binary ? "wb" : "w");
+  FILE *fp = Fopen(name.c_str(), binary ? "wb" : "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_UNV.cpp b/Geo/GModelIO_UNV.cpp
index f54504a57a7e7eb0aa93788fa98e05d3f231c9d3..2bb70db709f5e1106a4e06eb7a1ee9b90d328466 100644
--- a/Geo/GModelIO_UNV.cpp
+++ b/Geo/GModelIO_UNV.cpp
@@ -6,6 +6,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include "GModel.h"
+#include "OS.h"
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
@@ -16,7 +17,7 @@
 
 int GModel::readUNV(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -192,7 +193,7 @@ static std::string physicalName(GModel *m, int dim, int num)
 int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNodes,
                      double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_VRML.cpp b/Geo/GModelIO_VRML.cpp
index 99a7c9716e22976f0c1131e140cb380676aad341..8b354c5868941989a49a7936e7a28b275a0f67c6 100644
--- a/Geo/GModelIO_VRML.cpp
+++ b/Geo/GModelIO_VRML.cpp
@@ -6,6 +6,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include "GModel.h"
+#include "OS.h"
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
@@ -118,7 +119,7 @@ static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int r
 
 int GModel::readVRML(const std::string &name)
 {
-  FILE *fp = fopen(name.c_str(), "r");
+  FILE *fp = Fopen(name.c_str(), "r");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -189,7 +190,7 @@ int GModel::readVRML(const std::string &name)
 
 int GModel::writeVRML(const std::string &name, bool saveAll, double scalingFactor)
 {
-  FILE *fp = fopen(name.c_str(), "w");
+  FILE *fp = Fopen(name.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GModelIO_VTK.cpp b/Geo/GModelIO_VTK.cpp
index d4f1c096aa63b9242d6250aaa8c32f144c25119a..6fc8d12d118fea7b42991b1607fadf825b36beb2 100644
--- a/Geo/GModelIO_VTK.cpp
+++ b/Geo/GModelIO_VTK.cpp
@@ -4,6 +4,7 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "OS.h"
 #include "MPoint.h"
 #include "MLine.h"
 #include "MTriangle.h"
@@ -17,7 +18,7 @@
 int GModel::writeVTK(const std::string &name, bool binary, bool saveAll,
                      double scalingFactor, bool bigEndian)
 {
-  FILE *fp = fopen(name.c_str(), binary ? "wb" : "w");
+  FILE *fp = Fopen(name.c_str(), binary ? "wb" : "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
@@ -98,7 +99,7 @@ int GModel::writeVTK(const std::string &name, bool binary, bool saveAll,
 
 int GModel::readVTK(const std::string &name, bool bigEndian)
 {
-  FILE *fp = fopen(name.c_str(), "rb");
+  FILE *fp = Fopen(name.c_str(), "rb");
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index f9c424244be8c4eedeb629df5e23030b92b6e899..d1ad7e50e3a9fc9fe33269583ebcc6d32ac053d8 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -5,10 +5,11 @@
 //
 // Contributed by Cecile Piret
 
-#include "GmshConfig.h"
-#include "GRbf.h"
 #include <math.h>
 #include <vector>
+#include "GmshConfig.h"
+#include "GRbf.h"
+#include "OS.h"
 #include "fullMatrix.h"
 #include "Octree.h"
 #include "SPoint3.h"
@@ -59,19 +60,19 @@ static int SphereInEle(void *a, double*c)
 
 static void printNodes(std::set<MVertex *> myNodes)
 {
-  FILE * xyz = fopen("myNodes.pos","w");
+  FILE * xyz = Fopen("myNodes.pos","w");
   fprintf(xyz,"View \"\"{\n");
- for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end(); ++itv){
-   MVertex *v = *itv;
-   fprintf(xyz,"SP(%g,%g,%g){%d};\n", v->x(), v->y(), v->z(), v->getNum());
- }
- fprintf(xyz,"};\n");
- fclose(xyz);
+  for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end(); ++itv){
+    MVertex *v = *itv;
+    fprintf(xyz,"SP(%g,%g,%g){%d};\n", v->x(), v->y(), v->z(), v->getNum());
+  }
+  fprintf(xyz,"};\n");
+  fclose(xyz);
 }
 
 static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes)
 {
-  FILE *f = fopen ("UV.pos", "w");
+  FILE *f = Fopen ("UV.pos", "w");
   fprintf(f,"View  \" uv \" {\n");
 
   Msg::Info("*** RBF exporting 'UV.pos' ");
@@ -238,8 +239,8 @@ void GRbf::buildOctree(double radius)
 void GRbf::curvatureRBF(const fullMatrix<double> &cntrs,
 			fullMatrix<double> &curvature)
 {
-	
-	
+
+
   fullMatrix<double> extX, surf, sx,sy,sz, sxx,syy,szz, sxy,sxz,syz,sLap;
   setup_level_set(cntrs,normals,extX, surf);
 
@@ -253,9 +254,9 @@ void GRbf::curvatureRBF(const fullMatrix<double> &cntrs,
     double curv = -sLap(i,0)/norm_grad_s;
     curvature(i,0) = 0.5*fabs(curv)/sBox;
   }
-	
-	
-	
+
+
+
 }
 
 void GRbf::computeCurvature(const fullMatrix<double> &cntrs,
@@ -267,7 +268,7 @@ void GRbf::computeCurvature(const fullMatrix<double> &cntrs,
   evalRbfDer(1,extX,extX,surf,sx);
   evalRbfDer(2,extX,extX,surf,sy);
   evalRbfDer(3,extX,extX,surf,sz);
-  
+
   for (int i = 0; i < extX.size1(); i++) {
     double norm_grad_s = sqrt(sx(i,0)*sx(i,0)+sy(i,0)*sy(i,0)+sz(i,0)*sz(i,0));
     sx(i,0) = sx(i,0)/norm_grad_s;
@@ -296,7 +297,7 @@ void GRbf::computeCurvature(const fullMatrix<double> &cntrs,
 void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
 				 std::map<MVertex*, double> &rbf_curv)
 {
-  fullMatrix<double> extX, surf;  
+  fullMatrix<double> extX, surf;
   int numNodes = cntrs.size1();
   int numExtNodes = 3*numNodes;
   setup_level_set(cntrs,normals,extX, surf);
@@ -309,7 +310,7 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
     std::vector<int> &pts = nodesInSphere[i];
     int cluster_size = pts.size();
     fullMatrix<double> nodes_in_sph(3*cluster_size,3);
-    
+
     cluster_center(0,0) = extX(i,0);
     cluster_center(0,1) = extX(i,1);
     cluster_center(0,2) = extX(i,2);
@@ -345,11 +346,11 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
 	Dx(i+j*numNodes,pts[k]) = tempX(j,k);
 	Dy(i+j*numNodes,pts[k]) = tempY(j,k);
 	Dz(i+j*numNodes,pts[k]) = tempZ(j,k);
-	
+
 	Dx(i+j*numNodes,pts[k]+numNodes) = tempX(j,k+cluster_size);
 	Dy(i+j*numNodes,pts[k]+numNodes) = tempY(j,k+cluster_size);
 	Dz(i+j*numNodes,pts[k]+numNodes) = tempZ(j,k+cluster_size);
-	
+
 	Dx(i+j*numNodes,pts[k]+2*numNodes) = tempX(j,k+2*cluster_size);
 	Dy(i+j*numNodes,pts[k]+2*numNodes) = tempY(j,k+2*cluster_size);
 	Dz(i+j*numNodes,pts[k]+2*numNodes) = tempZ(j,k+2*cluster_size);
@@ -369,7 +370,7 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
     }
     sx.gemm(Dx,sx, 1.0, 0.0);
     sy.gemm(Dy,sy, 1.0, 0.0);
-    sz.gemm(Dz,sz, 1.0, 0.0); 
+    sz.gemm(Dz,sz, 1.0, 0.0);
 
     printf("sBox = %g ",sBox);
     for (int i = 0; i < numNodes; i++) {
@@ -506,7 +507,7 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
     cntrsPlus(i,0) = cntrs(i,0);
     cntrsPlus(i,1) = cntrs(i,1);
     cntrsPlus(i,2) = cntrs(i,2);
-  }	
+  }
   //Specifies the additional node position and its function value
   ONES(numNodes,0) = 1.0;
   cntrsPlus(numNodes,0) = cntrs(0,0)+10.*sBox;
@@ -514,7 +515,7 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
   cntrsPlus(numNodes,2) = cntrs(0,2)+10.*sBox;
 
   //printf("%g,%g,%g,%g;\n",sBox, cntrs(0,0), cntrs(0,1), cntrs(0,2));
-	
+
   evalRbfDer(1,cntrsPlus,cntrs,ONES,sx);
   evalRbfDer(2,cntrsPlus,cntrs,ONES,sy);
   evalRbfDer(3,cntrsPlus,cntrs,ONES,sz);
@@ -590,23 +591,23 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
 void GRbf::RbfLapSurface_global_projection2(const fullMatrix<double> &cntrs,
 										const fullMatrix<double> &normals,
 										fullMatrix<double> &Oper) {
-	
+
 	int numNodes = cntrs.size1();
 	int nnTot = 3*numNodes;
 	Oper.resize(numNodes,numNodes);
-	
+
 	fullMatrix<double> Dx(numNodes,numNodes),Dy(numNodes,numNodes),Dz(numNodes,numNodes),
     PDx(numNodes,numNodes),PDy(numNodes,numNodes),PDz(numNodes,numNodes),
     PDxx(numNodes,numNodes),PDyy(numNodes,numNodes),PDzz(numNodes,numNodes);
-	
-	fullMatrix<double> sx(nnTot,1),sy(nnTot,1),sz(nnTot,1); 
+
+	fullMatrix<double> sx(nnTot,1),sy(nnTot,1),sz(nnTot,1);
 	fullMatrix<double> extX(nnTot,3), surf(nnTot,1);
-	
+
 	//Stage 1 : The Arbitrary surface
 	setup_level_set(cntrs,normals,extX,surf);
 	if (!isLocal) extendedX = extX;
 	if (!isLocal) surfInterp = surf;
-	
+
 	//Computes the normal vectors to the surface at each node
 	evalRbfDer(1,extX,extX,surf,sx);
 	evalRbfDer(2,extX,extX,surf,sy);
@@ -621,7 +622,7 @@ void GRbf::RbfLapSurface_global_projection2(const fullMatrix<double> &cntrs,
 	RbfOp(1,cntrs,cntrs,Dx);
 	RbfOp(2,cntrs,cntrs,Dy);
 	RbfOp(3,cntrs,cntrs,Dz);
-	
+
 	// Fills up the operator matrix
 	for (int i=0;i<numNodes ; ++i){
 		for (int j=0;j<numNodes ; ++j){
@@ -638,8 +639,8 @@ void GRbf::RbfLapSurface_global_projection2(const fullMatrix<double> &cntrs,
 			Oper(i,j) = PDxx(i,j)+PDyy(i,j)+PDzz(i,j);
 		}
 	}
-	
-	
+
+
 }
 
 
@@ -723,7 +724,7 @@ void GRbf::RbfLapSurface_local_CPM(bool isLow,
     LocalOper.setAll(0.0);
     if (isLow) RbfLapSurface_global_CPM_low(nodes_in_sph,local_normals,LocalOper);
     else       RbfLapSurface_global_CPM_high_2(nodes_in_sph,local_normals,LocalOper);
-	  
+
     for (int j = 0; j < nbp; j++){
       Oper(i,pts[j])=LocalOper(0,j);
       Oper(i,pts[j]+numNodes)=LocalOper(0,j+nbp);
diff --git a/Geo/GeoStringInterface.cpp b/Geo/GeoStringInterface.cpp
index b171589ecb5b9efb3757a89c3ff117ef02281779..e449f04569139970945f8aeb6c69db0f90e55e6a 100644
--- a/Geo/GeoStringInterface.cpp
+++ b/Geo/GeoStringInterface.cpp
@@ -51,7 +51,7 @@ void add_infile(std::string text, std::string fileName, bool forceDestroy)
               return;
           }
         }
-        FILE *fp = fopen(newFileName.c_str(), "w");
+        FILE *fp = Fopen(newFileName.c_str(), "w");
         if(!fp) {
           Msg::Error("Unable to open file '%s'", newFileName.c_str());
           return;
@@ -69,14 +69,14 @@ void add_infile(std::string text, std::string fileName, bool forceDestroy)
 #if defined(HAVE_PARSER)
   std::string tmpFileName = CTX::instance()->homeDir + CTX::instance()->tmpFileName;
   FILE *gmsh_yyin_old = gmsh_yyin;
-  if(!(gmsh_yyin = fopen(tmpFileName.c_str(), "w"))) {
+  if(!(gmsh_yyin = Fopen(tmpFileName.c_str(), "w"))) {
     Msg::Error("Unable to open temporary file '%s'", tmpFileName.c_str());
     gmsh_yyin = gmsh_yyin_old;
     return;
   }
   fprintf(gmsh_yyin, "%s\n", text.c_str());
   fclose(gmsh_yyin);
-  gmsh_yyin = fopen(tmpFileName.c_str(), "r");
+  gmsh_yyin = Fopen(tmpFileName.c_str(), "r");
   while(!feof(gmsh_yyin)) {
     gmsh_yyparse();
   }
@@ -91,7 +91,7 @@ void add_infile(std::string text, std::string fileName, bool forceDestroy)
   GModel::current()->importGEOInternals();
   CTX::instance()->mesh.changed = ENT_ALL;
 
-  FILE *fp = fopen(fileName.c_str(), "a");
+  FILE *fp = Fopen(fileName.c_str(), "a");
   if(!fp) {
     Msg::Error("Unable to open file '%s'", fileName.c_str());
     return;
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index 93943e26168426f63e106714d42c958b7497dd62..20e19a48f2fa6218564006083ec8e70ca8cd6d7d 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -12,6 +12,7 @@
 #include "MEdge.h"
 #include "boundaryLayersData.h"
 #include "Field.h"
+#include "OS.h"
 
 SVector3 interiorNormal (SPoint2 p1, SPoint2 p2, SPoint2 p3)
 {
@@ -256,7 +257,7 @@ static void treat3Connections (GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2
   for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
 	 _columns->_normals.lower_bound(e3);
        itm != _columns->_normals.upper_bound(e3); ++itm) N3.push_back(itm->second);
-  
+
   SVector3 x1,x2;
   if (N1.size() == 2){
   }
@@ -369,7 +370,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
       SVector3 v02 = interiorNormal (p0,p2,p1);
       _columns->_normals.insert(std::make_pair(me02,v02));
     }
-    
+
     MEdge me21(v2,v1);
     if (allEdges.find(me21) != allEdges.end()){
       SVector3 v21 = interiorNormal (p2,p1,p0);
@@ -388,7 +389,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
     std::vector<SVector3> _dirs;
     // get all vertices that are connected  to that
     // vertex among all boundary layer vertices !
-    
+
     for (std::multimap<MVertex*,MVertex*>::iterator itm =
            _columns->_non_manifold_edges.lower_bound(*it);
          itm != _columns->_non_manifold_edges.upper_bound(*it); ++itm)
@@ -496,7 +497,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
       else {
 	MVertex *current = *it;
 	reparamMeshVertexOnFace(current,gf,p);
-	
+
 	int nbCol = 100;
 	std::vector<MVertex*> _column;
 	std::vector<SMetric3> _metrics;
@@ -505,7 +506,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
 	SPoint3 _close;
 	//double _current_distance = 0.;
 	while(1){
-	  
+
 	  SMetric3 m;
 	  double metric[3];
 	  double l;
@@ -547,7 +548,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
   }
   // DEBUG STUFF
 
-  FILE *f = fopen ("test.pos","w");
+  FILE *f = Fopen ("test.pos","w");
   fprintf(f,"View \"\" {\n");
   for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
     MVertex *v = *it;
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 53e6c0b404fb05ee578365b94fbbd809fcce534d..16f8fe1b2870599685b3fc4cded22ac8a8b55c82 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -13,6 +13,7 @@
 #include "gmshLevelset.h"
 #include "fullMatrix.h"
 #include "GModel.h"
+#include "OS.h"
 #include "MElement.h"
 #include "Numeric.h"
 #include "cartesian.h"
@@ -223,7 +224,7 @@ inline double evalRadialFnDer(int p, int index, double dx, double dy, double dz,
 
 inline void printNodes(fullMatrix<double> &myNodes, fullMatrix<double> &surf)
 {
-  FILE * xyz = fopen("myNodes.pos","w");
+  FILE * xyz = Fopen("myNodes.pos","w");
   fprintf(xyz,"View \"\"{\n");
   for(int itv = 1; itv != myNodes.size1(); ++itv){
     fprintf(xyz,"SP(%g,%g,%g){%g};\n", myNodes(itv,0), myNodes(itv,1),
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 72a54accbc001bd8934bb3a8d55c7cd5fab9a54c..e5c01e67a3e4a687eae97959d9dc9bccb49409bf 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -7,6 +7,7 @@
 #include <math.h>
 #include <stdio.h>
 #include "GmshMessage.h"
+#include "OS.h"
 #include "robustPredicates.h"
 #include "Numeric.h"
 #include "BDS.h"
@@ -16,7 +17,7 @@
 
 void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace *gf)
 {
-  FILE *f = fopen(iii, "w");
+  FILE *f = Fopen(iii, "w");
   if(!f) return;
   fprintf(f, "View \"scalar\" {\n");
   std::list<BDS_Face*>::iterator tit = t.begin();
@@ -533,7 +534,7 @@ void recur_tag(BDS_Face *t, BDS_GeomEntity *g)
 	//	recur_tag(t->e2->otherFace(t), g);
       }
       if(!t->e3->g && t->e3->numfaces() == 2) {
-	_stack.push(t->e3->otherFace(t));	
+	_stack.push(t->e3->otherFace(t));
 	//	recur_tag(t->e3->otherFace(t), g);
       }
     }
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index b1d8f0e09757f8eead36ac8fc39a33bafdfd701e..43ed51d7e5e71c0dc6562cf25b9677f8b5cedb7f 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -13,6 +13,7 @@
 #include "GFace.h"
 #include "GFaceCompound.h"
 #include "GModel.h"
+#include "OS.h"
 #include "Field.h"
 #include "MElement.h"
 #include "MElementOctree.h"
@@ -719,7 +720,7 @@ void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
 	SVector3 t2 (v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
 	t1.normalize();
 	t2.normalize();
-	double _angle = angle (t1,t2);	
+	double _angle = angle (t1,t2);
 	//        double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
         crossField2d::normalizeAngle (_angle);
         for (int i=0;i<2;i++){
@@ -766,7 +767,7 @@ void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
 inline double myAngle (const SVector3 &a, const SVector3 &b, const SVector3 &d){
   double cosTheta = dot(a,b);
   double sinTheta = dot(crossprod(a,b),d);
-  return atan2 (sinTheta,cosTheta);  
+  return atan2 (sinTheta,cosTheta);
 }
 
 
@@ -790,10 +791,10 @@ void backgroundMesh::propagateCrossField(GFace *_gf)
         reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
 	Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
 	SVector3 t1 = der.first();
-	SVector3 t2 = der.second(); 
+	SVector3 t2 = der.second();
 	SVector3 n = crossprod(t1,t2);
 	n.normalize();
-	SVector3 d1(v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z()); 
+	SVector3 d1(v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
 	t1.normalize();
 	d1.normalize();
 	double _angle = myAngle (t1,d1,n);
@@ -980,7 +981,7 @@ double backgroundMesh::getAngle(double u, double v, double w) const
 void backgroundMesh::print(const std::string &filename, GFace *gf,
                            const std::map<MVertex*,double> &_whatToPrint) const
 {
-  FILE *f = fopen (filename.c_str(),"w");
+  FILE *f = Fopen (filename.c_str(),"w");
   fprintf(f,"View \"Background Mesh\"{\n");
   for(unsigned int i=0;i<_triangles.size();i++){
     MVertex *v1 = _triangles[i]->getVertex(0);
@@ -996,7 +997,7 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
               v3->x(),v3->y(),v3->z(),itv1->second,itv2->second,itv3->second);
     }
     else {
-      
+
       GPoint p1 = gf->point(SPoint2(v1->x(),v1->y()));
       GPoint p2 = gf->point(SPoint2(v2->x(),v2->y()));
       GPoint p3 = gf->point(SPoint2(v3->x(),v3->y()));
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index c7b584cc61ce179df80f99e92b430763b6243e78..152e5db5e862f46696d9c2ef9338e2d0af3e5400 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -578,7 +578,7 @@ void Centerline::computeRadii()
 
 void Centerline::buildKdTree()
 {
-  FILE * f = fopen("myPOINTS.pos","w");
+  FILE * f = Fopen("myPOINTS.pos","w");
   fprintf(f, "View \"\"{\n");
 
   int nbPL = 3;  //10 points per line
@@ -1033,7 +1033,7 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad)
       // // theCut.insert(newCut.begin(),newCut.end());
       // char name[256];
       // sprintf(name, "myCUT-%d.pos", step);
-      // FILE * f2 = fopen(name,"w");
+      // FILE * f2 = Fopen(name,"w");
       // fprintf(f2, "View \"\"{\n");
       // std::set<MEdge,Less_Edge>::iterator itp =  newCut.begin();
       // while (itp != newCut.end()){
@@ -1282,7 +1282,7 @@ void Centerline::computeCrossField(double x,double y,double z,
 
 void Centerline::printSplit() const
 {
-  FILE * f = fopen("mySPLIT.pos","w");
+  FILE * f = Fopen("mySPLIT.pos","w");
   fprintf(f, "View \"\"{\n");
   for(unsigned int i = 0; i < edges.size(); ++i){
     std::vector<MLine*> lines = edges[i].lines;
@@ -1297,7 +1297,7 @@ void Centerline::printSplit() const
   fprintf(f,"};\n");
   fclose(f);
 
-  // FILE * f3 = fopen("myJUNCTIONS.pos","w");
+  // FILE * f3 = Fopen("myJUNCTIONS.pos","w");
   // fprintf(f3, "View \"\"{\n");
   //  std::set<MVertex*>::const_iterator itj = junctions.begin();
   //  while (itj != junctions.end()){
@@ -1310,7 +1310,7 @@ void Centerline::printSplit() const
   // fprintf(f3,"};\n");
   // fclose(f3);
 
-  FILE * f4 = fopen("myRADII.pos","w");
+  FILE * f4 = Fopen("myRADII.pos","w");
   fprintf(f4, "View \"\"{\n");
   for(unsigned int i = 0; i < lines.size(); ++i){
     MLine *l = lines[i];
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index aa42ddeacd6da401e4ef943f6f6f4caea0801a55..183da42e308f63ecb092614783e0d151181a85e6 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -20,7 +20,7 @@
 #include "Numeric.h"
 #include "robustPredicates.h"
 #include "BackgroundMesh.h"
-
+#include "OS.h"
 #include "GPoint.h"
 #include "GFace.h"
 #include "GEdgeCompound.h"
@@ -569,7 +569,7 @@ void DocRecord::voronoiCell(PointNumero pt, std::vector<SPoint2> &pts) const
 
 void DocRecord::makePosView(std::string fileName, GFace *gf)
 {
-  FILE *f = fopen(fileName.c_str(),"w");
+  FILE *f = Fopen(fileName.c_str(),"w");
    if (_adjacencies){
     fprintf(f,"View \"voronoi\" {\n");
     for(PointNumero i = 0; i < numPoints; i++) {
@@ -605,7 +605,7 @@ void DocRecord::makePosView(std::string fileName, GFace *gf)
 
 void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf, GEdge *ge)
 {
-  FILE *f = fopen(fileName.c_str(),"w");
+  FILE *f = Fopen(fileName.c_str(),"w");
    if (_adjacencies){
     fprintf(f,"View \"medial axis\" {\n");
     for(PointNumero i = 0; i < numPoints; i++) {
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index d87557af7d2519d56b0ced690a4e6b722deb578e..45bc5e09e7773aa68ebe959e887c766a9fac4adf 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -401,9 +401,9 @@ static void PrintMesh2dStatistics(GModel *m)
 
   FILE *statreport = 0;
   if(CTX::instance()->createAppendMeshStatReport == 1)
-    statreport = fopen(CTX::instance()->meshStatReportFileName.c_str(), "w");
+    statreport = Fopen(CTX::instance()->meshStatReportFileName.c_str(), "w");
   else if(CTX::instance()->createAppendMeshStatReport == 2)
-    statreport = fopen(CTX::instance()->meshStatReportFileName.c_str(), "a");
+    statreport = Fopen(CTX::instance()->meshStatReportFileName.c_str(), "a");
   else return;
 
   double worst = 1, best = 0, avg = 0;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 59f54a2bb052d762b817af7be708a26cbae48e35..71dde47beab12c4ec0c400c37f38fc48854ef893 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1183,7 +1183,7 @@ void printJacobians(GModel *m, const char *nm)
   const int n = 100;
   double D[n][n], X[n][n], Y[n][n], Z[n][n];
 
-  FILE *f = fopen(nm,"w");
+  FILE *f = Fopen(nm,"w");
   fprintf(f,"View \"\"{\n");
   for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
     for(unsigned int j = 0; j < (*it)->triangles.size(); j++){
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index 88f46c70cb7087823073322f838792c442f54235..bba0a7b1719a9c862c43c2be78b60112357119ae 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -12,6 +12,7 @@
 #include "meshGFaceDelaunayInsertion.h"
 #include "MTetrahedron.h"
 #include "directions3D.h"
+#include "OS.h"
 
 #if defined(HAVE_PETSC)
 #include "dofManager.h"
@@ -986,7 +987,7 @@ void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){
   std::map<MVertex*, STensor3>::iterator iter = crossField.find(beginV);
   STensor3 bCross = iter->second;
 
-  FILE *fi = fopen ("cross_recur.pos","w");
+  FILE *fi = Fopen ("cross_recur.pos","w");
   fprintf(fi,"View \"\"{\n");
   fprintf(fi,"SP(%g,%g,%g) {%g};\n",beginV->x(),beginV->y(),beginV->z(), 0.0);
   int count = 0;
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index b4be724d61516e358c935902c9c21a9aabd0ee92..303ea7b9c452f2a639966cc7093bf095b5646b5c 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -11,6 +11,7 @@
 
 #if defined(HAVE_SOLVER)
 
+#include "OS.h"
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
@@ -49,7 +50,7 @@ static double shapeMeasure(MElement *e)
   return d1;
 }
 
-void highOrderSmoother::moveTo(MVertex *v,  
+void highOrderSmoother::moveTo(MVertex *v,
                                    const std::map<MVertex*, SVector3> &m) const
 {
   std::map<MVertex*,SVector3>::const_iterator it = m.find(v);
@@ -58,7 +59,7 @@ void highOrderSmoother::moveTo(MVertex *v,
     v->y() = it->second.y();
     v->z() = it->second.z();
   }
-} 
+}
 
 void highOrderSmoother::moveToStraightSidedLocation(MVertex *v) const
 {
@@ -82,7 +83,7 @@ void highOrderSmoother::moveToTargetLocation(MElement *e) const
     moveToTargetLocation(e->getVertex(i));
 }
 
-void highOrderSmoother::updateTargetLocation(MVertex*v, const SPoint3 &p3, 
+void highOrderSmoother::updateTargetLocation(MVertex*v, const SPoint3 &p3,
                                                  const SPoint2 &p2)
 {
   v->x() = p3.x();
@@ -99,7 +100,7 @@ struct p2data{
   MVertex *n12;
   fullMatrix<double> *m1, *m2;
   highOrderSmoother *s;
-  p2data(GFace *_gf, MTriangle *_t1, MTriangle *_t2, MVertex *_n12, 
+  p2data(GFace *_gf, MTriangle *_t1, MTriangle *_t2, MVertex *_n12,
          highOrderSmoother *_s)
     : gf(_gf), t1(_t1), t2(_t2), n12(_n12), s(_s)
   {
@@ -109,7 +110,7 @@ struct p2data{
     m1 = new fullMatrix<double>(3 * t1->getNumVertices(),
                                 3 * t1->getNumVertices());
     m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
-                                3 * t2->getNumVertices()); 
+                                3 * t2->getNumVertices());
     SElement se1(t1), se2(t2);
     el.elementMatrix(&se1, *m1);
     el.elementMatrix(&se2, *m2);
@@ -139,7 +140,7 @@ struct pNdata{
     m1 = new  fullMatrix<double>(3 * t1->getNumVertices(),
                                  3 * t1->getNumVertices());
     m2 = new  fullMatrix<double>(3 * t2->getNumVertices(),
-                                 3 * t2->getNumVertices()); 
+                                 3 * t2->getNumVertices());
     SElement se1(t1), se2(t2);
     el.elementMatrix(&se1, *m1);
     el.elementMatrix(&se2, *m2);
@@ -153,7 +154,7 @@ struct pNdata{
   }
 };
 
-static double _DeformationEnergy(MElement *e, 
+static double _DeformationEnergy(MElement *e,
                                  fullMatrix<double> *K,
                                  highOrderSmoother *s)
 {
@@ -183,13 +184,13 @@ static double _DeformationEnergy(MElement *e,
 
 static double _DeformationEnergy(p2data *p)
 {
-  return _DeformationEnergy(p->t1, p->m1, p->s) + 
+  return _DeformationEnergy(p->t1, p->m1, p->s) +
     _DeformationEnergy(p->t2, p->m2, p->s);
 }
 
 static double _DeformationEnergy(pNdata *p)
 {
-  return _DeformationEnergy(p->t1, p->m1, p->s) + 
+  return _DeformationEnergy(p->t1, p->m1, p->s) +
     _DeformationEnergy(p->t2, p->m2, p->s);
 }
 
@@ -207,7 +208,7 @@ static double _function_p2tB(fullVector<double> &x, void *data)
   double E = _DeformationEnergy(p);
   p->n12->x() = xx;
   p->n12->y() = yy;
-  p->n12->z() = zz;  
+  p->n12->z() = zz;
   return E;
 }
 
@@ -221,10 +222,10 @@ static double _function_p2t(fullVector<double> &x, void *data)
   p->n12->x() = gp12.x();
   p->n12->y() = gp12.y();
   p->n12->z() = gp12.z();
-  double q = std::min(shapeMeasure(p->t1), shapeMeasure(p->t2));      
+  double q = std::min(shapeMeasure(p->t1), shapeMeasure(p->t2));
   p->n12->x() = xx;
   p->n12->y() = yy;
-  p->n12->z() = zz;  
+  p->n12->z() = zz;
   return -q;
 }
 
@@ -243,12 +244,12 @@ static double _function_pNt(fullVector<double> &x, void *data)
     p->n[i]->y() = gp12.y();
     p->n[i]->z() = gp12.z();
   }
-  double q = std::min(shapeMeasure(p->t1), shapeMeasure(p->t2));      
+  double q = std::min(shapeMeasure(p->t1), shapeMeasure(p->t2));
   for (unsigned int i = 0; i < p->n.size(); i++){
     p->n[i]->x() = xx[i];
     p->n[i]->y() = yy[i];
     p->n[i]->z() = zz[i];
-  }  
+  }
   return -q;
 }
 
@@ -272,7 +273,7 @@ static double _function_pNtB(fullVector<double> &x, void *data)
     p->n[i]->x() = xx[i];
     p->n[i]->y() = yy[i];
     p->n[i]->z() = zz[i];
-  }  
+  }
   return E;
 }
 
@@ -291,7 +292,7 @@ void getDistordedElements(const std::vector<MElement*> &v,
   }
 }
 
-void addOneLayer(const std::vector<MElement*> &v, 
+void addOneLayer(const std::vector<MElement*> &v,
                  std::vector<MElement*> &d,
                  std::vector<MElement*> &layer)
 {
@@ -327,7 +328,7 @@ void addOneLayer(const std::vector<MElement*> &v,
 
 void _printJacobiansAtNodes (const char * name, std::vector<MElement*> &v)
 {
-  FILE *fd = fopen (name,"w");
+  FILE *fd = Fopen (name,"w");
   fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$ElementNodeData\n1\n\"det J\"\n1\n0.0\n3\n1\n1\n%d\n", (int) v.size());
   for (int i=0; i<v.size(); i++){
     const polynomialBasis* pb = v[i]->getFunctionSpace();
@@ -338,9 +339,9 @@ void _printJacobiansAtNodes (const char * name, std::vector<MElement*> &v)
       const double _w = pb->points(j,2);
       double jac[3][3];
       double J = v[i]->getJacobian (_u,_v,_w, jac);
-      fprintf(fd," %g",J);    
+      fprintf(fd," %g",J);
     }
-    fprintf(fd,"\n");    
+    fprintf(fd,"\n");
   }
   fprintf(fd,"$EndElementNodeData\n");
   fclose(fd);
@@ -384,14 +385,14 @@ void highOrderSmoother::smooth(GRegion *gr)
 
 
 
-void highOrderSmoother::optimize(GFace * gf, 
+void highOrderSmoother::optimize(GFace * gf,
 				 edgeContainer &edgeVertices,
 				 faceContainer &faceVertices)
 {
   //  if (gf->geomType() != GEntity::Plane) return;
 
     std::vector<MElement*> bad;
-    
+
 
   while (1) {
     // relocate the vertices using elliptic smoother
@@ -400,17 +401,17 @@ void highOrderSmoother::optimize(GFace * gf,
     //      findOptimalLocationsPN(gf,this);
     //      findOptimalLocationsPN(gf,this);
     //    }
-    // then try to swap for better configurations  
+    // then try to swap for better configurations
 
     //smooth(gf, !CTX::instance()->mesh.highOrderNoMetric);
-    
-    
+
+
     //    for (int i=0;i<100;i++){
-      //      int nbSwap = 
+      //      int nbSwap =
 	//	swapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
 	//      printf("%d swaps\n",nbSwap);
     //    }
-    
+
      smooth(gf,true);
     // smooth(gf,true);
     // smooth(gf,true);
@@ -429,8 +430,8 @@ void highOrderSmoother::optimize(GFace * gf,
   }
 }
 
-void highOrderSmoother::computeMetricVector(GFace *gf, 
-                                                MElement *e, 
+void highOrderSmoother::computeMetricVector(GFace *gf,
+                                                MElement *e,
                                                 fullMatrix<double> &J,
                                                 fullMatrix<double> &JT,
                                                 fullVector<double> &D)
@@ -438,8 +439,8 @@ void highOrderSmoother::computeMetricVector(GFace *gf,
   int nbNodes = e->getNumVertices();
   for (int j = 0; j < nbNodes; j++){
     SPoint2 param;
-    reparamMeshVertexOnFace(e->getVertex(j), gf, param);  
-    Pair<SVector3,SVector3> der = gf->firstDer(param);    
+    reparamMeshVertexOnFace(e->getVertex(j), gf, param);
+    Pair<SVector3,SVector3> der = gf->firstDer(param);
     int XJ = j;
     int YJ = j + nbNodes;
     int ZJ = j + 2 * nbNodes;
@@ -451,7 +452,7 @@ void highOrderSmoother::computeMetricVector(GFace *gf,
     J(XJ,VJ) = der.second().x();
     J(YJ,VJ) = der.second().y();
     J(ZJ,VJ) = der.second().z();
-    
+
     JT(UJ,XJ) = der.first().x();
     JT(UJ,YJ) = der.first().y();
     JT(UJ,ZJ) = der.first().z();
@@ -459,12 +460,12 @@ void highOrderSmoother::computeMetricVector(GFace *gf,
     JT(VJ,YJ) = der.second().y();
     JT(VJ,ZJ) = der.second().z();
 
-    GPoint gp = gf->point(param);    
+    GPoint gp = gf->point(param);
     SVector3 ss = getSSL(e->getVertex(j));
     D(XJ) = gp.x() - ss.x();
     D(YJ) = gp.y() - ss.y();
     D(ZJ) = gp.z() - ss.z();
-  } 
+  }
 }
 
 void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
@@ -479,10 +480,10 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
 #endif
   dofManager<double> myAssembler(lsys);
   elasticityTerm El(0, 1.0, .33, getTag());
-  
+
   std::vector<MElement*> layer, v;
   double minD;
-  
+
   getDistordedElements(all, 0.99, v, minD);
 
   const int nbLayers = 3; //2, originally :)
@@ -519,10 +520,10 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
       // printf("%d %d %d v\n", i, j, v[i]->getNumVertices());
       its = _straightSidedLocation.find(vert);
       if (its == _straightSidedLocation.end()){
-        _straightSidedLocation[vert] = 
-          SVector3(vert->x(), vert->y(), vert->z());     
-        _targetLocation[vert] = 
-          SVector3(vert->x(), vert->y(), vert->z());     
+        _straightSidedLocation[vert] =
+          SVector3(vert->x(), vert->y(), vert->z());
+        _targetLocation[vert] =
+          SVector3(vert->x(), vert->y(), vert->z());
       }
       else{
         vert->x() = its->second.x();
@@ -535,7 +536,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
       }
     }
   }
-  
+
   // number the other DOFs
   for (unsigned int i = 0; i < v.size(); i++){
     for (int j = 0; j < v[i]->getNumVertices(); j++){
@@ -543,9 +544,9 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
       myAssembler.numberVertex(vert, 0, getTag());
       myAssembler.numberVertex(vert, 1, getTag());
       verticesToMove.insert(vert);
-    } 
+    }
   }
-  
+
   double dx0 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
   double dx = dx0;
   Msg::Debug(" dx0 = %12.5E", dx0);
@@ -560,8 +561,8 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
 
   for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
     SPoint2 param;
-    reparamMeshVertexOnFace(*it, gf, param);  
-    GPoint gp = gf->point(param);    
+    reparamMeshVertexOnFace(*it, gf, param);
+    GPoint gp = gf->point(param);
     if ((*it)->onWhat()->dim() == 2){
       (*it)->x() = gp.x();
       (*it)->y() = gp.y();
@@ -572,14 +573,14 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
       SVector3 p =  _targetLocation[(*it)];
       (*it)->x() = p.x();
       (*it)->y() = p.y();
-      (*it)->z() = p.z();      
+      (*it)->z() = p.z();
     }
   }
   delete lsys;
 }
 
-double highOrderSmoother::smooth_metric_(std::vector<MElement*>  & v, 
-                                         GFace *gf, 
+double highOrderSmoother::smooth_metric_(std::vector<MElement*>  & v,
+                                         GFace *gf,
                                          dofManager<double> &myAssembler,
                                          std::set<MVertex*> &verticesToMove,
                                          elasticityTerm &El)
@@ -624,11 +625,11 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*>  & v,
     myAssembler.systemSolve();
     // for all element, compute detJ at integration points --> material law
     // end while convergence -------------------------------------------------------
-  
+
     for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
       if ((*it)->onWhat()->dim() == 2){
 	SPoint2 param;
-	reparamMeshVertexOnFace((*it), gf, param);  
+	reparamMeshVertexOnFace((*it), gf, param);
 	SPoint2 dparam;
 	myAssembler.getDofValue((*it), 0, getTag(), dparam[0]);
 	myAssembler.getDofValue((*it), 1, getTag(), dparam[1]);
@@ -637,10 +638,10 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*>  & v,
 	(*it)->setParameter(0, newp.x());
 	(*it)->setParameter(1, newp.y());
       }
-    }    
+    }
     myAssembler.systemClear();
   }
-  
+
   return dx;
 }
 
@@ -654,7 +655,7 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
 {
 #ifdef HAVE_PETSC
   // assume that the mesh is OK, yet already curved
-  linearSystemPETSc<double> *lsys = new  linearSystemPETSc<double>;    
+  linearSystemPETSc<double> *lsys = new  linearSystemPETSc<double>;
   dofManager<double> myAssembler(lsys);
   elasticityMixedTerm El_mixed (0, 1.0, .333, getTag());
   elasticityTerm El (0, 1.0, .333, getTag());
@@ -662,8 +663,8 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
   std::set<MVertex*> _vertices;
 
   //+++++++++ Boundary Conditions & Numbering +++++++++++++++++++++++++++++++
-  // fix all dof that correspond to vertices on the boundary 
-  // the value is equal 
+  // fix all dof that correspond to vertices on the boundary
+  // the value is equal
   for (unsigned int i = 0; i < v.size(); i++){
     for (int j = 0; j < v[i]->getNumVertices(); j++){
       MVertex *vert = v[i]->getVertex(j);
@@ -699,7 +700,7 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
       }
     }
     if (_dim == 2)myAssembler.fixVertex(vert, 2, getTag(), 0);
-    // number vertices 
+    // number vertices
     myAssembler.numberVertex(vert, 0, getTag());
     myAssembler.numberVertex(vert, 1, getTag());
     myAssembler.numberVertex(vert, 2, getTag());
@@ -720,15 +721,15 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
     // solve the system
     lsys->systemSolve();
   }
-  
+
   //+++++++++ Move vertices @ maximum ++++++++++++++++++++++++++++++++++++
-  FILE *fd = fopen ("d.msh","w");
+  FILE *fd = Fopen ("d.msh","w");
   fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$NodeData\n1\n\"tr(sigma)\"\n1\n0.0\n3\n1\n3\n%d\n", (int) _vertices.size());
   for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
     double ax, ay, az;
     myAssembler.getDofValue(*it, 0, getTag(), ax);
     myAssembler.getDofValue(*it, 1, getTag(), ay);
-    myAssembler.getDofValue(*it, 2, getTag(), az);    
+    myAssembler.getDofValue(*it, 2, getTag(), az);
     (*it)->x() += max_incr*ax;
     (*it)->y() += max_incr*ay;
     (*it)->z() += max_incr*az;
@@ -744,7 +745,7 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
   /*
     std::map<MVertex*,double> facts;
     while(1){
-    int count = 0; 
+    int count = 0;
     for (unsigned int i = 0; i < v.size(); i++){
     double disto = v[i]->distoShapeMeasure();
     if (disto < thres){
@@ -758,7 +759,7 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
     double ax, ay, az;
     myAssembler.getDofValue(vert, 0, getTag(), ax);
     myAssembler.getDofValue(vert, 1, getTag(), ay);
-    myAssembler.getDofValue(vert, 2, getTag(), az);    
+    myAssembler.getDofValue(vert, 2, getTag(), az);
     vert->x() -= .1*ax;
     vert->y() -= .1*ay;
     vert->z() -= .1*az;
@@ -774,14 +775,14 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
   while(1){
     std::vector<MElement*> disto;
     double minD;
-    getDistordedElements(v, 0.5, disto, minD);    
+    getDistordedElements(v, 0.5, disto, minD);
     if (minD < thres){
       percentage -= 10.;
       for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
 	double ax, ay, az;
 	myAssembler.getDofValue(*it, 0, getTag(), ax);
 	myAssembler.getDofValue(*it, 1, getTag(), ay);
-	myAssembler.getDofValue(*it, 2, getTag(), az);    
+	myAssembler.getDofValue(*it, 2, getTag(), az);
 	(*it)->x() -= .1*ax;
 	(*it)->y() -= .1*ay;
 	(*it)->z() -= .1*az;
@@ -789,14 +790,14 @@ double highOrderSmoother::apply_incremental_displacement (double max_incr,
     }
     else break;
   }
-    
+
   delete lsys;
   return percentage;
 #endif
   return 0.0;
 }
 
-int highOrderSmoother::smooth_with_mixed_formulation (std::vector<MElement*> &all, 
+int highOrderSmoother::smooth_with_mixed_formulation (std::vector<MElement*> &all,
 						      double alpha)
 {
   int ITER = 0;
@@ -818,18 +819,18 @@ int highOrderSmoother::smooth_with_mixed_formulation (std::vector<MElement*> &al
     else if (percentage_of_what_is_left == 100.)break;
   }
 
-  getDistordedElements(all, 0.3, disto, minD);    
+  getDistordedElements(all, 0.3, disto, minD);
   Msg::Info("iter %d : %d elements / %d distorted  min Disto = %g ",ITER,
-	    all.size(), disto.size(), minD);       
+	    all.size(), disto.size(), minD);
 
   if (0 && percentage < 0.99){
     char NN[256];
     sprintf(NN,"smoothing-%d.msh",ITER++);
-    double percentage_of_what_is_left = apply_incremental_displacement (1.0,all, true, alpha,NN,disto);    
+    double percentage_of_what_is_left = apply_incremental_displacement (1.0,all, true, alpha,NN,disto);
     percentage += (1.-percentage) * percentage_of_what_is_left/100.;
     Msg::Info("The mixed smoother was able to do %3d percent of the motion",(int)(percentage*100.));
   }
-  
+
   return 1;
 }
 
@@ -839,7 +840,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
 #if defined(HAVE_TAUCS)
   linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
 #elif defined(HAVE_PETSC)
-  linearSystemPETSc<double> *lsys = new  linearSystemPETSc<double>;    
+  linearSystemPETSc<double> *lsys = new  linearSystemPETSc<double>;
 #else
   linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
   lsys->setNoisy(1);
@@ -848,7 +849,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
 #endif
   dofManager<double> myAssembler(lsys);
   elasticityTerm El(0, 1.0, .333, getTag());
-  
+
   std::vector<MElement*> layer, v;
   double minD;
 
@@ -909,7 +910,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
 
   for (unsigned int i = 0; i < v.size(); i++)
     moveToStraightSidedLocation(v[i]);
-  
+
   // number the other DOFs
   for (unsigned int i = 0; i < v.size(); i++){
     for (int j = 0; j < v[i]->getNumVertices(); j++){
@@ -919,7 +920,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
       myAssembler.numberVertex(vert, 2, getTag());
       // gather all vertices that are supposed to move
       verticesToMove[vert] = SVector3(0.0, 0.0, 0.0);
-    } 
+    }
   }
 
   //  Msg::Info("%d vertices FIXED %d NUMBERED", myAssembler.sizeOfF()
@@ -936,7 +937,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
     // solve the system
     lsys->systemSolve();
   }
-  
+
   for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
     double ax, ay, az;
     myAssembler.getDofValue(it->first, 0, getTag(), ax);
@@ -948,14 +949,14 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
   }
 
   // delete matrices and vectors
-  
+
   delete lsys;
 }
 
 void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
                                       std::vector<MElement*>& old_elems,
                                       GFace* gf) {
- 
+
   //  printf("Smoothing a cavity...\n");
   //  printf("Old elems : %d and %d\n", old_elems[0]->getNum(), old_elems[1]->getNum());
   //  printf("Cavity elems : %d and %d\n", cavity[0]->getNum(), cavity[1]->getNum());
@@ -996,7 +997,7 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
       }
     }
   }
-  
+
   //printf("%d vertices \n", _displ.size());
 
   std::set<MVertex*>::iterator it;
@@ -1014,9 +1015,9 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
       its = _straightSidedLocation.find(vert);
       if (its == _straightSidedLocation.end()) {
 	//printf("SETTING LOCATIONS for %d\n",vert->getNum());
-        _straightSidedLocation[vert] = 
-          SVector3(vert->x(), vert->y(), vert->z());     
-        _targetLocation[vert] = 
+        _straightSidedLocation[vert] =
+          SVector3(vert->x(), vert->y(), vert->z());
+        _targetLocation[vert] =
           SVector3(vert->x(), vert->y(), vert->z());
       }else {
         vert->x() = its->second.x();
@@ -1039,12 +1040,12 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
       myAssembler.numberVertex(vert, 0, getTag());
       myAssembler.numberVertex(vert, 1, getTag());
       verticesToMove.insert(vert);
-    } 
+    }
   }
 
   //Msg::Info("%d vertices FIXED %d NUMBERED\n", myAssembler.sizeOfF()
   //          , myAssembler.sizeOfR());
-  
+
   double dx0 = smooth_metric_(cavity, gf, myAssembler, verticesToMove, El);
   double dx = dx0;
   //printf(" dx0 = %12.5E\n", dx0);
@@ -1059,8 +1060,8 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
 
   for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
     SPoint2 param;
-    reparamMeshVertexOnFace(*it, gf, param);  
-    GPoint gp = gf->point(param);    
+    reparamMeshVertexOnFace(*it, gf, param);
+    GPoint gp = gf->point(param);
     if ((*it)->onWhat()->dim() == 2){
       (*it)->x() = gp.x();
       (*it)->y() = gp.y();
@@ -1071,7 +1072,7 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
       SVector3 p =  _targetLocation[(*it)];
       (*it)->x() = p.x();
       (*it)->y() = p.y();
-      (*it)->z() = p.z();      
+      (*it)->z() = p.z();
     }
     //printf("  Moving %d to %g %g %g\n", (*it)->getNum(),_targetLocation[(*it)][0], _targetLocation[(*it)][1],_targetLocation[(*it)][2] );
   }
@@ -1095,8 +1096,8 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
 */
 
 /*
-  n34 is the new vertex, we'd like to put it at the 
-  best place i.e. at a place where it maximizes the 
+  n34 is the new vertex, we'd like to put it at the
+  best place i.e. at a place where it maximizes the
   minimal smoothness of the neighboring elements.
 
   One can actually maximize this value or one can use
@@ -1110,8 +1111,8 @@ static void getNodesP2(const MEdge &me, MTriangle *t1, MTriangle *t2,
                        MVertex * &n13){
 
   n1 = me.getVertex(0);
-  n2 = me.getVertex(1);    
-  
+  n2 = me.getVertex(1);
+
   if (t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
   else if (t1->getVertex(1) != n1 && t1->getVertex(1) != n2) n3 = t1->getVertex(1);
   else if (t1->getVertex(2) != n1 && t1->getVertex(2) != n2) n3 = t1->getVertex(2);
@@ -1129,7 +1130,7 @@ static void getNodesP2(const MEdge &me, MTriangle *t1, MTriangle *t2,
     n23 = t1->getVertex((iStart+1)%3 + 3);
     n12 = t1->getVertex((iStart+2)%3 + 3);
   }
-  
+
   if (t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
   else if (t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
   else if (t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
@@ -1156,8 +1157,8 @@ static void getNodesPN (const MEdge &me, MTriangle *t1, MTriangle *t2,
                         std::vector<MVertex *> &n13){
 
   n1 = me.getVertex(0);
-  n2 = me.getVertex(1);    
-  
+  n2 = me.getVertex(1);
+
   int order = t1->getPolynomialOrder();
 
   if (t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
@@ -1183,7 +1184,7 @@ static void getNodesPN (const MEdge &me, MTriangle *t1, MTriangle *t2,
     int start12 = 3 + ((iStart+2)%3) * (order-1);
     for (int i=order-2;i>=0;i--)n12.push_back(t1->getVertex(start12+i));
   }
-  
+
   if (t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
   else if (t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
   else if (t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
@@ -1206,7 +1207,7 @@ static void getNodesPN (const MEdge &me, MTriangle *t1, MTriangle *t2,
   }
 }
 
-static double surfaceTriangleUV(SPoint2 &v1, SPoint2 &v2, SPoint2 &v3)           
+static double surfaceTriangleUV(SPoint2 &v1, SPoint2 &v2, SPoint2 &v3)
 {
   const double v12[2] = {v2.x() - v1.x(),v2.y() - v1.y()};
   const double v13[2] = {v3.x() - v1.x(),v3.y() - v1.y()};
@@ -1227,7 +1228,7 @@ struct swap_triangles_p2
   double evalConfiguration (GFace *gf, SPoint2 & p) const
   {
     GPoint gp34 = gf->point(p);
-    MFaceVertex _test (gp34.x(),gp34.y(),gp34.z(),gf,p.x(),p.y());        
+    MFaceVertex _test (gp34.x(),gp34.y(),gp34.z(),gf,p.x(),p.y());
     std::vector<MVertex *>vv;
     vv.push_back(n13);vv.push_back(&_test);vv.push_back(n14);
     MTriangleN t3_test (n1,n3,n4,vv,2,t1->getNum(),t1->getPartition());
@@ -1239,15 +1240,15 @@ struct swap_triangles_p2
     return std::min(q3,q4);
   }
 
-  MVertex *optimalLocation (GFace *gf, 
+  MVertex *optimalLocation (GFace *gf,
                             SPoint2 &p3,
                             SPoint2 &p4) const
   {
     SPoint2 p34_linear = (p3+p4)*.5;
-    
+
     GPoint gp34 = gf->point(p34_linear);
     MFaceVertex *_test = new MFaceVertex (gp34.x(),gp34.y(),gp34.z(),
-                                          gf,p34_linear.x(),p34_linear.y());        
+                                          gf,p34_linear.x(),p34_linear.y());
     std::vector<MVertex *>vv;
     vv.push_back(n13);vv.push_back(_test);vv.push_back(n14);
     MTriangleN t3_test (n1,n3,n4,vv,2,t1->getNum(),t1->getPartition());
@@ -1312,7 +1313,7 @@ struct swap_triangles_p2
            n2->getNum(),n3->getNum(),n4->getNum(),
            n12->getNum(),n23->getNum(),n13->getNum(),n24->getNum(),n14->getNum());
   }
-  
+
 };
 
 struct swap_triangles_pN
@@ -1335,7 +1336,7 @@ struct swap_triangles_pN
     std::pair<MVertex*, MVertex*> _temp2(_temp.getMinVertex(), _temp.getMaxVertex());
     edgeVertices.erase(_temp2);
     // do the same for internal vertices !!
-  } 
+  }
 
   swap_triangles_pN(const MEdge &me, MTriangle *_t1, MTriangle *_t2, GFace *gf,
                     edgeContainer &_edgeVertices,
@@ -1372,7 +1373,7 @@ struct swap_triangles_pN
     old_elems.push_back(t1);
     old_elems.push_back(t2);
     //s->smooth_cavity(cavity, old_elems, gf);
-      
+
     const double qnew1 = shapeMeasure(t3);
     const double qnew2 = shapeMeasure(t4);
 
@@ -1389,15 +1390,15 @@ struct swap_triangles_pN
   }
 };
 
-static int optimalLocationP2_(GFace *gf, 
+static int optimalLocationP2_(GFace *gf,
                               const MEdge &me,
-                              MTriangle *t1, MTriangle *t2, 
+                              MTriangle *t1, MTriangle *t2,
                               highOrderSmoother *s)
 {
   double qini = std::min(shapeMeasure(t1),shapeMeasure(t2));
 
   if (qini > 0.6) return 0;
-  
+
   MVertex *n1,*n2,*n3=0,*n4=0,*n12,*n14,*n24,*n23,*n13;
   getNodesP2 (me,t1,t2,n1,n2,n3,n4,n12,n14,n24,n23,n13);
   SPoint2 p1,p2,p3,p4,p12;
@@ -1423,8 +1424,8 @@ static int optimalLocationP2_(GFace *gf,
   double u1 = gf->parBounds(0).high();
   double v0 = gf->parBounds(1).low();
   double v1 = gf->parBounds(1).high();
-  if (pp(0) < u0 || pp(0) > u1 || pp(1) < v0 || pp(1) > v1)return 0;   
-  
+  if (pp(0) < u0 || pp(0) > u1 || pp(1) < v0 || pp(1) > v1)return 0;
+
   GPoint gp12 = gf->point(SPoint2(pp(0),pp(1)));
   n12->x() = gp12.x();
   n12->y() = gp12.y();
@@ -1510,7 +1511,7 @@ static int findOptimalLocationsPN(GFace *gf,highOrderSmoother *s)
   e2t_cont adj;
   buildEdgeToTriangle(gf->triangles, adj);
   int N=0;
-  
+
   for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
     if (it->second.second)
       N += optimalLocationPN_(gf,it->first, dynamic_cast<MTriangle*>(it->second.first),
@@ -1519,7 +1520,7 @@ static int findOptimalLocationsPN(GFace *gf,highOrderSmoother *s)
   return N;
 }
 
-static int swapHighOrderTriangles(GFace *gf, 
+static int swapHighOrderTriangles(GFace *gf,
                                   edgeContainer &edgeVertices,
                                   faceContainer &faceVertices,
                                   highOrderSmoother *s)
@@ -1613,7 +1614,7 @@ static int swapHighOrderTriangles(GFace *gf,
           v_removed.insert(*vit);
 	  //delete  *vit;
       }
-      
+
       swap_triangles_pN &sw = (swap_triangles_pN &) *itp;
       sw.cleanupDeletedEdge (false);
 
@@ -1626,7 +1627,7 @@ static int swapHighOrderTriangles(GFace *gf,
   for (unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
     if (v_removed.find(gf->mesh_vertices[i]) == v_removed.end()){
       mesh_vertices2.push_back(gf->mesh_vertices[i]);
-    } 
+    }
   }
 
   gf->mesh_vertices.clear();
@@ -1638,7 +1639,7 @@ static int swapHighOrderTriangles(GFace *gf,
     }
     else {
       delete gf->triangles[i];
-    }    
+    }
   }
   gf->triangles.clear();
   gf->triangles = triangles2;
@@ -1700,14 +1701,14 @@ static int swapHighOrderTriangles(GFace *gf)
     }
     ++itp;
   }
-    
+
   for (unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
     if (v_removed.find(gf->mesh_vertices[i]) == v_removed.end()){
       mesh_vertices2.push_back(gf->mesh_vertices[i]);
     }
     else {
       delete gf->mesh_vertices[i];
-    }    
+    }
   }
   gf->mesh_vertices = mesh_vertices2;
 
@@ -1717,14 +1718,14 @@ static int swapHighOrderTriangles(GFace *gf)
     }
     else {
       delete gf->triangles[i];
-    }    
+    }
   }
   //  printf("replacing %d by %d\n",gf->triangles.size(),triangles2.size());
   gf->triangles = triangles2;
   return nbSwap;
 }
 
-void  highOrderSmoother::swap(GFace *gf, 
+void  highOrderSmoother::swap(GFace *gf,
                                   edgeContainer &edgeVertices,
                                   faceContainer &faceVertices)
 {
@@ -1758,10 +1759,10 @@ highOrderSmoother::highOrderSmoother (GModel *gm)
 {
   _clean();
   // compute straigh sided positions that are actually X,Y,Z positions
-  // that are NOT always on curves and surfaces 
+  // that are NOT always on curves and surfaces
 
   // compute stright sided positions of vertices that are classified on model edges
-  for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){    
+  for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
     for (int i=0;i<(*it)->lines.size();i++){
       MLine *l = (*it)->lines[i];
       int N = l->getNumVertices()-2;
@@ -1790,8 +1791,8 @@ highOrderSmoother::highOrderSmoother (GModel *gm)
 	  const double t1 = fs->points(j, 0);
 	  const double t2 = fs->points(j, 1);
 	  SPoint3 pc = face.interpolate(t1, t2);
-	  _straightSidedLocation [t->getVertex(j)] = 
-	    SVector3(pc.x(),pc.y(),pc.z()); 
+	  _straightSidedLocation [t->getVertex(j)] =
+	    SVector3(pc.x(),pc.y(),pc.z());
 	}
       }
     }
@@ -1804,8 +1805,8 @@ highOrderSmoother::highOrderSmoother (GModel *gm)
 	  const double t1 = fs->points(j, 0);
 	  const double t2 = fs->points(j, 1);
 	  SPoint3 pc = face.interpolate(t1, t2);
-	  _straightSidedLocation [q->getVertex(j)] = 
-	    SVector3(pc.x(),pc.y(),pc.z()); 
+	  _straightSidedLocation [q->getVertex(j)] =
+	    SVector3(pc.x(),pc.y(),pc.z());
 	}
       }
     }
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index f65f5edda0b060b4afdc182c8ed3b034b2641b48..49200f113ee893b0d95b89a5dd1bb6a54a252e50 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -644,7 +644,7 @@ double highOrderTools::apply_incremental_displacement (double max_incr,
   }
 
   //+++++++++ Move vertices @ maximum ++++++++++++++++++++++++++++++++++++
-  FILE *fd = fopen ("d.msh","w");
+  FILE *fd = Fopen ("d.msh","w");
   fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$NodeData\n1\n\"tr(sigma)\"\n1\n0.0\n3\n1\n3\n%d\n", (int) _vertices.size());
   for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
     double ax, ay, az;
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 5555a1afaccade5a0da6765907ea82a8d38f3b37..b790eb7ce6cde30ee80ea8749b42e98ec81d9135 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -15,6 +15,7 @@
 #include "Context.h"
 #include "STensor3.h"
 #include "Field.h"
+#include "OS.h"
 
 #define SQU(a)      ((a)*(a))
 
@@ -105,7 +106,7 @@ static double F_Lc_aniso(GEdge *ge, double t)
 #else
   bool blf = false;
 #endif
-  
+
 
   GPoint p = ge->point(t);
   SMetric3 lc_here;
@@ -125,7 +126,7 @@ static double F_Lc_aniso(GEdge *ge, double t)
   if (blf && !blf->isEdgeBL(ge->tag())){
     SMetric3 lc_bgm;
     blf->computeFor1dMesh ( p.x(), p.y(), p.z() , lc_bgm );
-    lc_here = intersection_conserveM1 (lc_here, lc_bgm );			    
+    lc_here = intersection_conserveM1 (lc_here, lc_bgm );
   }
 #endif
 
@@ -301,7 +302,7 @@ static void printFandPrimitive(int tag, std::vector<IntPoint> &Points)
 {
   char name[256];
   sprintf(name, "line%d.dat", tag);
-  FILE *f = fopen(name, "w");
+  FILE *f = Fopen(name, "w");
   if(!f) return;
   double l = 0;
   for (unsigned int i = 0; i < Points.size(); i++){
@@ -493,7 +494,7 @@ void meshGEdge::operator() (GEdge *ge)
 
 #if defined(HAVE_ANN)
   if (blf && !blf->isEdgeBL(ge->tag()))
-    {      
+    {
       GVertex *g0 = ge->getBeginVertex();
       GVertex *g1 = ge->getEndVertex();
       BoundaryLayerColumns* _columns = ge->model()->getColumns();
@@ -513,7 +514,7 @@ void meshGEdge::operator() (GEdge *ge)
       }
       if (blf->isVertexBL(g1->tag())){
 	  _columns->addColumn(t*-1.0, v1,invert,_metrics);
-      }      
+      }
     }
 #endif
   ge->meshStatistics.status = GEdge::DONE;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 25579c3a62f9e8df01ddac4406d0d3d405c2c1d6..797402350799164274a94894667dfe6008421fe4 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -645,7 +645,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
   std::list<GEdge*> embedded_edges = gf->embeddedEdges();
   edges.insert(edges.begin(), embedded_edges.begin(),embedded_edges.end());
   std::list<GEdge*>::iterator ite = edges.begin();
-  FILE *ff2 = fopen ("tato.pos","w");
+  FILE *ff2 = Fopen ("tato.pos","w");
   fprintf(ff2,"View \" \"{\n");
   std::set<MVertex*> verts;
   while(ite != edges.end()){
@@ -762,7 +762,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
   std::list<GEdge*> hop;
   std::set<MEdge,Less_Edge>::iterator it =  bedges.begin();
 
-  FILE *ff = fopen ("toto.pos","w");
+  FILE *ff = Fopen ("toto.pos","w");
   fprintf(ff,"View \" \"{\n");
   for (; it != bedges.end(); ++it){
     ne.lines.push_back(new MLine (it->getVertex(0),it->getVertex(1)));
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 3c67143c9889f1492a051d5927841bc8bfa72027..aca7e4f7ff78011888a6043332c98a4d50b09897 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -45,7 +45,7 @@ static bool isBoundary(MTri3 *t, double limit_, int &active)
 template <class ITERATOR>
 void _printTris(char *name, ITERATOR it,  ITERATOR end, bidimMeshData & data, bool param=true)
 {
-  FILE *ff = fopen (name,"w");
+  FILE *ff = Fopen (name,"w");
   fprintf(ff,"View\"test\"{\n");
   while ( it != end ){
     MTri3 *worst = *it;
@@ -211,9 +211,9 @@ void circumCenterMetric(MTriangle *base, const double *metric, bidimMeshData & d
                         double *x, double &Radius2)
 {
   // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
-  int index0 = data.getIndex (base->getVertex(0)); 
-  int index1 = data.getIndex (base->getVertex(1)); 
-  int index2 = data.getIndex (base->getVertex(2)); 
+  int index0 = data.getIndex (base->getVertex(0));
+  int index1 = data.getIndex (base->getVertex(1));
+  int index2 = data.getIndex (base->getVertex(2));
   double pa[2] = {data.Us[index0], data.Vs[index0] };
   double pb[2] = {data.Us[index1], data.Vs[index1] };
   double pc[2] = {data.Us[index2], data.Vs[index2] };
@@ -283,9 +283,9 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base,
   double x[2], Radius2;
   double metric[3];
   if (!metricb){
-    int index0 = data.getIndex (base->getVertex(0)); 
-    int index1 = data.getIndex (base->getVertex(1)); 
-    int index2 = data.getIndex (base->getVertex(2)); 
+    int index0 = data.getIndex (base->getVertex(0));
+    int index1 = data.getIndex (base->getVertex(1));
+    int index2 = data.getIndex (base->getVertex(2));
     double pa[2] = {(data.Us[index0] +data.Us[index1] + data.Us[index2]) / 3.,
   		  (data.Vs[index0] +data.Vs[index1] + data.Vs[index2]) / 3.};
     buildMetric(gf, pa, metric);
@@ -330,9 +330,9 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, bidimMeshData * data, GF
     }
     else {
 
-      int index0 = data->getIndex (base->getVertex(0)); 
-      int index1 = data->getIndex (base->getVertex(1)); 
-      int index2 = data->getIndex (base->getVertex(2)); 
+      int index0 = data->getIndex (base->getVertex(0));
+      int index1 = data->getIndex (base->getVertex(1));
+      int index2 = data->getIndex (base->getVertex(2));
       double p1[2] = {data->Us[index0], data->Vs[index0] };
       double p2[2] = {data->Us[index1], data->Vs[index1] };
       double p3[2] = {data->Us[index2], data->Vs[index2] };
@@ -386,9 +386,9 @@ int MTri3::inCircumCircle(const double *p) const
 int inCircumCircle(MTriangle *base,
                    const double *p, const double *param , bidimMeshData & data)
 {
-  int index0 = data.getIndex (base->getVertex(0)); 
-  int index1 = data.getIndex (base->getVertex(1)); 
-  int index2 = data.getIndex (base->getVertex(2)); 
+  int index0 = data.getIndex (base->getVertex(0));
+  int index1 = data.getIndex (base->getVertex(1));
+  int index2 = data.getIndex (base->getVertex(2));
   double pa[2] = {data.Us[index0], data.Vs[index0] };
   double pb[2] = {data.Us[index1], data.Vs[index1] };
   double pc[2] = {data.Us[index2], data.Vs[index2] };
@@ -495,7 +495,7 @@ void recurFindCavityAniso(GFace *gf,
     MTri3 *neigh = t->getNeigh(i) ;
     edgeXface exf (t, i);
     // take care of untouchable internal edges
-    std::set<MEdge,Less_Edge>::iterator it = data.internalEdges.find(MEdge(exf.v[0],exf.v[1])); 
+    std::set<MEdge,Less_Edge>::iterator it = data.internalEdges.find(MEdge(exf.v[0],exf.v[1]));
     if (!neigh || it != data.internalEdges.end())
       shell.push_back(exf);
     else  if (!neigh->isDeleted()){
@@ -511,9 +511,9 @@ void recurFindCavityAniso(GFace *gf,
 
 bool circUV(MTriangle *t, bidimMeshData & data, double *res, GFace *gf)
 {
-  int index0 = data.getIndex (t->getVertex(0)); 
-  int index1 = data.getIndex (t->getVertex(1)); 
-  int index2 = data.getIndex (t->getVertex(2)); 
+  int index0 = data.getIndex (t->getVertex(0));
+  int index1 = data.getIndex (t->getVertex(1));
+  int index2 = data.getIndex (t->getVertex(2));
   double u1[3] = {data.Us[index0], data.Vs[index0], 0 };
   double u2[3] = {data.Us[index1], data.Vs[index1], 0 };
   double u3[3] = {data.Us[index2], data.Vs[index2], 0 };
@@ -527,9 +527,9 @@ bool invMapUV(MTriangle *t, double *p, bidimMeshData & data,
   double mat[2][2];
   double b[2];
 
-  int index0 = data.getIndex (t->getVertex(0)); 
-  int index1 = data.getIndex (t->getVertex(1)); 
-  int index2 = data.getIndex (t->getVertex(2)); 
+  int index0 = data.getIndex (t->getVertex(0));
+  int index1 = data.getIndex (t->getVertex(1));
+  int index2 = data.getIndex (t->getVertex(2));
 
   double u0 = data.Us[index0];
   double v0 = data.Vs[index0];
@@ -559,9 +559,9 @@ bool invMapUV(MTriangle *t, double *p, bidimMeshData & data,
 
 inline double getSurfUV(MTriangle *t, bidimMeshData & data)
 {
-  int index0 = data.getIndex (t->getVertex(0)); 
-  int index1 = data.getIndex (t->getVertex(1)); 
-  int index2 = data.getIndex (t->getVertex(2)); 
+  int index0 = data.getIndex (t->getVertex(0));
+  int index1 = data.getIndex (t->getVertex(1));
+  int index2 = data.getIndex (t->getVertex(2));
 
   double u1 = data.Us[index0];
   double v1 = data.Vs[index0];
@@ -607,9 +607,9 @@ bool insertVertexB (std::list<edgeXface> &shell,
   bool onePointIsTooClose = false;
   while (it != shell.end()){
     MTriangle *t = new MTriangle(it->v[0], it->v[1], v);
-    int index0 = data.getIndex (t->getVertex(0)); 
-    int index1 = data.getIndex (t->getVertex(1)); 
-    int index2 = data.getIndex (t->getVertex(2)); 
+    int index0 = data.getIndex (t->getVertex(0));
+    int index1 = data.getIndex (t->getVertex(1));
+    int index2 = data.getIndex (t->getVertex(2));
     const double ONE_THIRD = 1./3.;
     double lc = ONE_THIRD * (data.vSizes[index0] +data.vSizes[index1] +data.vSizes[index2]);
     double lcBGM = ONE_THIRD * (data.vSizesBGM[index0] +data.vSizesBGM[index1] +data.vSizesBGM[index2]);
@@ -734,9 +734,9 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2], bidimMeshData & data,
   int ITER = 0;
   while (1){
     N_SEARCH ++ ;
-    int index0 = data.getIndex (t->tri()->getVertex(0)); 
-    int index1 = data.getIndex (t->tri()->getVertex(1)); 
-    int index2 = data.getIndex (t->tri()->getVertex(2)); 
+    int index0 = data.getIndex (t->tri()->getVertex(0));
+    int index1 = data.getIndex (t->tri()->getVertex(1));
+    int index2 = data.getIndex (t->tri()->getVertex(2));
     SPoint3 q2((data.Us[index0] +data.Us[index1] +data.Us[index2])/ 3.0,
 	       (data.Vs[index0] +data.Vs[index1] +data.Vs[index2])/ 3.0,
 	       0);
@@ -813,7 +813,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
     //      printf("cocuou\n");
     ptin = search4Triangle (worst, center, data, AllTris,uv, oneNewTriangle ? true : false);
       if (ptin) {
-	recurFindCavityAniso(gf, shell, cavity, metric, center, ptin, data);    
+	recurFindCavityAniso(gf, shell, cavity, metric, center, ptin, data);
       }
   }
 
@@ -832,11 +832,11 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
     int index2 = data.getIndex(ptin->tri()->getVertex(2));
     lc1 = (1. - uv[0] - uv[1]) * data.vSizes[index0] + uv[0] * data.vSizes[index1] + uv[1] * data.vSizes[index2];
     lc = BGM_MeshSize(gf, center[0], center[1], p.x(), p.y(), p.z());
-    
+
     //SMetric3 metr = BGM_MeshMetric(gf, center[0], center[1], p.x(), p.y(), p.z());
     //                               vMetricsBGM.push_back(metr);
-    data.addVertex ( v ,  center[0], center[1], lc1, lc ); 
-    
+    data.addVertex ( v ,  center[0], center[1], lc1, lc );
+
     //    double t1 = Cpu();
 
     if(!p.succeeded() || !insertVertexB
@@ -871,7 +871,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
 }
 
 void gmshRuppert(GFace *gf,  double minqual, int MAXPNT,
-		 std::map<MVertex* , MVertex*>* equivalence,  
+		 std::map<MVertex* , MVertex*>* equivalence,
 		 std::map<MVertex*, SPoint2> * parametricCoordinates){
   MTri3::radiusNorm =3;
 
@@ -920,7 +920,7 @@ void gmshRuppert(GFace *gf,  double minqual, int MAXPNT,
 
 
 void bowyerWatson(GFace *gf, int MAXPNT,
-		  std::map<MVertex* , MVertex*>* equivalence,  
+		  std::map<MVertex* , MVertex*>* equivalence,
 		  std::map<MVertex*, SPoint2> * parametricCoordinates)
 {
   std::set<MTri3*,compareTri3Ptr> AllTris;
@@ -993,7 +993,7 @@ void bowyerWatson(GFace *gf, int MAXPNT,
     }
   }
 #endif
-  transferDataStructure(gf, AllTris, DATA);  
+  transferDataStructure(gf, AllTris, DATA);
 }
 
 /*
@@ -1206,7 +1206,7 @@ void optimalPointFrontalB (GFace *gf,
 }
 
 void bowyerWatsonFrontal(GFace *gf,
-		  std::map<MVertex* , MVertex*>* equivalence,  
+		  std::map<MVertex* , MVertex*>* equivalence,
 		  std::map<MVertex*, SPoint2> * parametricCoordinates)
 {
   std::set<MTri3*,compareTri3Ptr> AllTris;
@@ -1381,7 +1381,7 @@ void optimalPointFrontalQuadB (GFace *gf,
 }
 
 void buildBackGroundMesh (GFace *gf,
-		  std::map<MVertex* , MVertex*>* equivalence,  
+		  std::map<MVertex* , MVertex*>* equivalence,
 		  std::map<MVertex*, SPoint2> * parametricCoordinates)
 {
   //  printf("build bak mesh\n");
@@ -1430,7 +1430,7 @@ void buildBackGroundMesh (GFace *gf,
 }
 
 void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
-		  std::map<MVertex* , MVertex*>* equivalence,  
+		  std::map<MVertex* , MVertex*>* equivalence,
 		  std::map<MVertex*, SPoint2> * parametricCoordinates)
 {
 
@@ -1568,7 +1568,7 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
 }
 
 void bowyerWatsonParallelograms(GFace *gf,
-		  std::map<MVertex* , MVertex*>* equivalence,  
+		  std::map<MVertex* , MVertex*>* equivalence,
 		  std::map<MVertex*, SPoint2> * parametricCoordinates)
 {
   std::set<MTri3*,compareTri3Ptr> AllTris;
diff --git a/Mesh/meshGFaceElliptic.cpp b/Mesh/meshGFaceElliptic.cpp
index 6571842881f239a4755358204666238fd880006e..b2affefce27cf68c65eb69f1dba2154e6d57ab79 100644
--- a/Mesh/meshGFaceElliptic.cpp
+++ b/Mesh/meshGFaceElliptic.cpp
@@ -33,14 +33,14 @@
 #define TRAN_QUAD(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
     (1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4)
 
-
-static void printQuads(GFace *gf, fullMatrix<SPoint2> uv, std::vector<MQuadrangle*> quads,  int iter){
+static void printQuads(GFace *gf, fullMatrix<SPoint2> uv,
+                       std::vector<MQuadrangle*> quads,  int iter){
 
   if(!CTX::instance()->mesh.saveAll) return;
 
   char name[234];
   sprintf(name,"quadUV_%d_%d.pos", gf->tag(), iter);
-  FILE *f = fopen(name,"w");
+  FILE *f = Fopen(name,"w");
   fprintf(f,"View \"%s\" {\n",name);
 
   for (int i = 1; i < uv.size1()-1; i++)
@@ -52,7 +52,7 @@ static void printQuads(GFace *gf, fullMatrix<SPoint2> uv, std::vector<MQuadrangl
 
   char name3[234];
   sprintf(name3,"quadXYZ_%d_%d.pos", gf->tag(), iter);
-  FILE *f3 = fopen(name3,"w");
+  FILE *f3 = Fopen(name3,"w");
   fprintf(f3,"View \"%s\" {\n",name3);
   for (unsigned int i = 0; i < quads.size(); i++){
     quads[i]->writePOS(f3,true,false,false,false,false,false);
@@ -86,7 +86,7 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M
 
   char name[234];
   sprintf(name,"paramGrid_%d.pos", gf->tag());
-  FILE *f = fopen(name,"w");
+  FILE *f = Fopen(name,"w");
   fprintf(f,"View \"%s\" {\n",name);
 
   // for (unsigned int i = 0; i < p1.size(); i++)
@@ -106,7 +106,7 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M
 
   char name2[234];
   sprintf(name2,"paramEdges_%d.pos", gf->tag());
-  FILE *f2 = fopen(name2,"w");
+  FILE *f2 = Fopen(name2,"w");
   fprintf(f2,"View \"%s\" {\n",name2);
 
   for (unsigned int i = 0; i < e01.size(); i++){
@@ -143,7 +143,7 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M
 
   char name3[234];
   sprintf(name3,"quadXYZ_%d.pos", gf->tag());
-  FILE *f3 = fopen(name3,"w");
+  FILE *f3 = Fopen(name3,"w");
   fprintf(f3,"View \"%s\" {\n",name3);
   for (unsigned int i = 0; i < quads.size(); i++){
     quads[i]->writePOS(f3,true,false,false,false,false,false);
@@ -155,8 +155,8 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M
 
 }
 
-static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv, 
-			      std::vector<std::vector<MVertex*> > &tab, 
+static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv,
+			      std::vector<std::vector<MVertex*> > &tab,
 			      std::vector<MQuadrangle*> &newq,  std::vector<MVertex*> &newv){
 
   newq.clear();
@@ -164,17 +164,17 @@ static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv,
 
   int M = uv.size1();
   int N = uv.size2();
- 
+
   for (int i = 1; i < M-1; i++){
     for (int j = 0; j < N-1; j++){
       GPoint gp = gf->point(uv(i,j));
-      if (!gp.succeeded()){ 
+      if (!gp.succeeded()){
   	printf("** QUADS FROM UV new vertex not created p=%g %g \n", uv(i,j).x(), uv(i,j).y()); //exit(1);
       }
       tab[i][j] = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
     }
   }
-  for (int i = 1; i < M-1; i++) tab[i][N-1] =  tab[i][0]; 
+  for (int i = 1; i < M-1; i++) tab[i][N-1] =  tab[i][0];
 
   //create vertices
   for (int i = 1; i < M-1; i++)
@@ -188,7 +188,7 @@ static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv,
       newq.push_back(qnew);
     }
   }
-  
+
   return;
 
 }
@@ -235,10 +235,10 @@ static std::vector<MVertex*> saturateEdgeHarmonic (GFace *gf, SPoint2 p1, SPoint
 }
 
 static void transfiniteSmoother(GFace* gf,
-				fullMatrix<SPoint2> &uv, 
+				fullMatrix<SPoint2> &uv,
 				std::vector<std::vector<MVertex*> > &tab,
 				std::vector<MQuadrangle*> &newq,
-				std::vector<MVertex*> &newv, 
+				std::vector<MVertex*> &newv,
 				bool isPeriodic=false){
 
    int M = uv.size1();
@@ -257,18 +257,18 @@ static void transfiniteSmoother(GFace* gf,
   			       SQU(uv(i,jp1).y() - uv(i,jm1).y())) ;
   	double gamma = 0.25 * (SQU(uv(i+1,j).x() - uv(i-1,j).x()) +
   			       SQU(uv(i+1,j).y() - uv(i-1,j).y()));
-  	double beta = 0.0625 * 
+  	double beta = 0.0625 *
   	  ((uv(i+1,j).x() - uv(i-1,j).x()) * (uv(i,jp1).x() - uv(i,jm1).x()) +
   	   (uv(i+1,j).y() - uv(i-1,j).y()) * (uv(i,jp1).y() - uv(i,jm1).y()));
-  	double uk = 0.5 * 
-  	  (alpha * (uv(i+1,j).x() + uv(i-1,j).x()) + 
+  	double uk = 0.5 *
+  	  (alpha * (uv(i+1,j).x() + uv(i-1,j).x()) +
   	   gamma * (uv(i,jp1).x() + uv(i,jm1).x()) -
-  	   2. * beta * (uv(i+1,jp1).x() - uv(i-1,jp1).x() - 
+  	   2. * beta * (uv(i+1,jp1).x() - uv(i-1,jp1).x() -
   			uv(i+1,jm1).x() + uv(i-1,jm1).x())) / (alpha + gamma);
-  	double vk = 0.5 * 
-  	  (alpha * (uv(i+1,j).y() + uv(i-1,j).y()) + 
+  	double vk = 0.5 *
+  	  (alpha * (uv(i+1,j).y() + uv(i-1,j).y()) +
   	   gamma * (uv(i,jp1).y()+ uv(i,jm1).y()) -
-  	   2. * beta * (uv(i+1,jp1).y() - uv(i-1,jp1).y() - 
+  	   2. * beta * (uv(i+1,jp1).y() - uv(i-1,jp1).y() -
   			uv(i+1,jm1).y() + uv(i-1,jm1).y())) / (alpha + gamma);
   	uvold(i,j) = uv(i,j);
   	norm += sqrt((uv(i,j).x()-uk)*(uv(i,j).x()-uk)+(uv(i,j).y()-vk)*(uv(i,j).y()-vk));
@@ -293,7 +293,7 @@ static void transfiniteSmoother(GFace* gf,
   //final print
   createQuadsFromUV(gf, uv, tab, newq, newv);
   printQuads(gf, uv, newq, numSmooth);
- 
+
 }
 
 //elliptic surface grid generator
@@ -302,7 +302,7 @@ static void ellipticSmoother( GFace* gf,
 			      fullMatrix<SPoint2> &uv,
 			      std::vector<std::vector<MVertex*> > &tab,
 			      std::vector<MQuadrangle*> &newq,
-			      std::vector<MVertex*> &newv, 
+			      std::vector<MVertex*> &newv,
 			      bool isPeriodic=false){
 
   printQuads(gf, uv, newq, 0);
@@ -333,8 +333,8 @@ static void ellipticSmoother( GFace* gf,
 	double g11_bar = dot(xu,xu);
 	double g12_bar = dot(xu,xv);
 	double g22_bar = dot(xv,xv);
-	double dxsi = 1.; 
-	double deta = 1.; 
+	double dxsi = 1.;
+	double deta = 1.;
 	double u_xsi = (uv(i,jp1).x()-uv(i,j).x())/dxsi;
 	double u_eta = (uv(i+1,j).x()-uv(i,j).x())/deta;
 	double v_xsi = (uv(i,jp1).y()-uv(i,j).y())/dxsi;
@@ -342,7 +342,7 @@ static void ellipticSmoother( GFace* gf,
 	double g11 = g11_bar*u_xsi*u_xsi+2.*g12_bar*u_xsi*v_xsi+g22_bar*v_xsi*v_xsi;
 	double g12 = g11_bar*u_xsi*u_eta+g12_bar*(u_xsi*v_eta+u_eta*v_xsi)+g22_bar*v_xsi*v_eta;
 	double g22 = g11_bar*u_eta*u_eta+2.*g12_bar*u_eta*v_eta+g22_bar*v_eta*v_eta;
-	double jac = u_xsi*v_eta-u_eta*v_xsi; //sqrt(g11*g22-g12*g12);	
+	double jac = u_xsi*v_eta-u_eta*v_xsi; //sqrt(g11*g22-g12*g12);
 	double jac_bar = sqrt(g11_bar*g22_bar-g12_bar*g12_bar);
 	double g11_bar_u = 2.*dot(xu,xuu);
 	double g11_bar_v = 2.*dot(xu,xuv);
@@ -374,12 +374,12 @@ static void ellipticSmoother( GFace* gf,
 	uvold(i,N-1) = uvold(i,0);
       }
     }
-    
+
     //under-relaxation
     // double omega = 0.8;
     // for (int i=0; i<M; i++){
     //   for (int j = 0; j<N; j++){
-    // 	uv(i,j) = SPoint2(omega*uv(i,j).x() + (1.-omega)*uvold(i,j).x(), 
+    // 	uv(i,j) = SPoint2(omega*uv(i,j).x() + (1.-omega)*uvold(i,j).x(),
     // 			  omega*uv(i,j).y() + (1.-omega)*uvold(i,j).y());
     //   }
     // }
@@ -424,15 +424,15 @@ static void createRegularGrid (GFace *gf,
 			       std::vector<MVertex*> &e41, std::vector<SPoint2> &pe41,int sign41,
 			       std::vector<MQuadrangle*> &q,
 			       std::vector<MVertex*> &newv,
-			       fullMatrix<SPoint2> &uv, 
+			       fullMatrix<SPoint2> &uv,
 			       std::vector<std::vector<MVertex*> > &tab)
 {
- 
+
   int M = e23.size();
   int N = e12.size();
 
   uv.resize(M+2,N+2);
- 
+
   //std::vector<std::vector<MVertex*> > tab(M+2);
   tab.resize(M+2);
   for(int i = 0; i <= M+1; i++) tab[i].resize(N+2);
@@ -441,7 +441,7 @@ static void createRegularGrid (GFace *gf,
   tab[0][N+1] = v2;    uv(0,N+1) = c2;
   tab[M+1][N+1] = v3;  uv(M+1,N+1) = c3;
   tab[M+1][0] = v4;    uv(M+1,0) = c4;
- 
+
   for (int i=0;i<N;i++){
     tab[0][i+1]   = sign12 > 0 ? e12 [i] : e12 [N-i-1];
     tab[M+1][i+1] = sign34 < 0 ? e34 [i] : e34 [N-i-1];
@@ -511,15 +511,15 @@ static void createRegularGrid (GFace *gf,
 static void createRegularGridPeriodic  (GFace *gf,int sign2,
 					MVertex *v0, SPoint2 &c0,
 					MVertex *v2, SPoint2 &c2,
-					std::vector<MVertex*> &e02, std::vector<SPoint2> &pe02, 
-					std::vector<MVertex*> &e00, std::vector<SPoint2> &pe00, 
-					std::vector<MVertex*> &e22, std::vector<SPoint2> &pe22, 
+					std::vector<MVertex*> &e02, std::vector<SPoint2> &pe02,
+					std::vector<MVertex*> &e00, std::vector<SPoint2> &pe00,
+					std::vector<MVertex*> &e22, std::vector<SPoint2> &pe22,
 					std::vector<MQuadrangle*> &q,
 					std::vector<MVertex*> &newv,
-					fullMatrix<SPoint2> &uv, 
+					fullMatrix<SPoint2> &uv,
 					std::vector<std::vector<MVertex*> > &tab)
 {
- 
+
   int M = e02.size();
   int N = e00.size();
 
@@ -527,7 +527,7 @@ static void createRegularGridPeriodic  (GFace *gf,int sign2,
 
   char name3[234];
   sprintf(name3,"quadParam_%d.pos", gf->tag());
-  FILE *f3 = fopen(name3,"w");
+  FILE *f3 = Fopen(name3,"w");
   fprintf(f3,"View \"%s\" {\n",name3);
 
   tab.resize(M+2);
@@ -539,8 +539,8 @@ static void createRegularGridPeriodic  (GFace *gf,int sign2,
   tab[M+1][N+1] = v2;  uv(M+1,N+1) = c2;
   tab[M+1][0] = v2;    uv(M+1,0) = c2;
   for (int i=0;i<N;i++){
-    tab[0][i+1]   =   e00 [i];  uv(0,i+1) = pe00 [i]  ; 
-    tab[M+1][i+1] = (sign2 > 0) ?    e22 [i]  : e22 [N-i-1]  ;  
+    tab[0][i+1]   =   e00 [i];  uv(0,i+1) = pe00 [i]  ;
+    tab[M+1][i+1] = (sign2 > 0) ?    e22 [i]  : e22 [N-i-1]  ;
     uv(M+1,i+1)  =  (sign2 > 0) ?    pe22 [i] : pe22 [N-i-1] ;
   }
   for (int i=0;i<M;i++){
@@ -550,12 +550,12 @@ static void createRegularGridPeriodic  (GFace *gf,int sign2,
 
   //interior mesh_vertices
   for (int i=0;i<N;i++){
-    SPoint2 p0 =  pe00[i]  ; 
+    SPoint2 p0 =  pe00[i]  ;
     SPoint2 p2 =  (sign2>0) ?    pe22[i] : pe22 [N-i-1] ;
     std::vector<SPoint2> pei;
     std::vector<MVertex*> ei = saturateEdgeRegular(gf,p0,p2,M+1,pei);//M+1
     for (int j=0;j<M;j++){
-      SPoint2 pij = pei[j]; 
+      SPoint2 pij = pei[j];
       GPoint gp = gf->point(pij);
       if (!gp.succeeded()) printf("** INITIAL GRID new vertex not created p=%g %g \n", pij.x(), pij.y());
       tab[j+1][i+1] = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
@@ -567,7 +567,7 @@ static void createRegularGridPeriodic  (GFace *gf,int sign2,
   for (int i=0;i<N+1;i++)
     for (int j=1;j<M+1;j++)
       newv.push_back(tab[j][i]);
-      
+
   // create quads
   for (int i=0;i<=N;i++){
     for (int j=0;j<=M;j++){
@@ -594,12 +594,12 @@ static void updateFaceQuads(GFace *gf, std::vector<MQuadrangle*> &quads, std::ve
   for(unsigned int i = 0; i < newv.size(); i++){
     gf->mesh_vertices.push_back(newv[i]);
   }
-  
+
 }
 
-static bool computeRingVertices(GFace *gf, Centerline *center, 
-				std::vector<MVertex*> &vert1, std::vector<MVertex*> &vert2, 
-				int &N, int &M, int &close_ind, int &sign2, 
+static bool computeRingVertices(GFace *gf, Centerline *center,
+				std::vector<MVertex*> &vert1, std::vector<MVertex*> &vert2,
+				int &N, int &M, int &close_ind, int &sign2,
 				double &arc, double &length){
 
 #if defined(HAVE_ANN)
@@ -662,10 +662,10 @@ static bool computeRingVertices(GFace *gf, Centerline *center,
 
   SPoint2 p1; reparamMeshVertexOnFace(vert1[0], gf, p1);
   SPoint2 p2; reparamMeshVertexOnFace(vert2[close_ind], gf, p2);
-  GPoint gp_i(vert1[0]->x(), vert1[0]->y(), vert1[0]->z()); 
+  GPoint gp_i(vert1[0]->x(), vert1[0]->y(), vert1[0]->z());
   SPoint2 p1b; reparamMeshVertexOnFace(vert1[N/2], gf, p1b);
   SPoint2 p2b; reparamMeshVertexOnFace(vert2[(close_ind+N/2)%N], gf, p2b);
-  GPoint gp_ib(vert1[N/2]->x(), vert1[N/2]->y(), vert1[N/2]->z()); 
+  GPoint gp_ib(vert1[N/2]->x(), vert1[N/2]->y(), vert1[N/2]->z());
   length = 0.0;
   double lengthb = 0.0;
   int nb = 100;
@@ -674,8 +674,8 @@ static bool computeRingVertices(GFace *gf, Centerline *center,
     GPoint gp_ii = gf->point(pii);
     SPoint2 piib(p1b.x() + (double)(i+1)/nb*(p2b.x()-p1b.x()), p1b.y() + (double)(i+1)/nb*(p2b.y()-p1b.y()));
     GPoint gp_iib = gf->point(piib);
-    length  += gp_i.distance(gp_ii); 
-    lengthb += gp_ib.distance(gp_iib); 
+    length  += gp_i.distance(gp_ii);
+    lengthb += gp_ib.distance(gp_iib);
     gp_i = gp_ii;
     gp_ib = gp_iib;
   }
@@ -683,7 +683,7 @@ static bool computeRingVertices(GFace *gf, Centerline *center,
   arc = 2*M_PI*rad;
   double lc =  arc/N;
   M = (int)(0.5*(length+lengthb)/lc) ;
-  
+
   delete kdtree;
   delete[]index;
   delete[]dist;
@@ -720,7 +720,7 @@ bool createRegularTwoCircleGridPeriodic (Centerline *center, GFace *gf)
   std::vector<MVertex*> e00,e22;//edges without first and last vertex
   for (int i=1;i<N;i++) e00.push_back(vert1[i]);
   for (int i=1;i<N;i++) e22.push_back(vert2[(close_ind+i)%N]);
- 
+
   std::vector<SPoint2> pe00, pe22;
   for (unsigned int i = 0; i < e00.size(); i++){
      SPoint2 p00; reparamMeshVertexOnFace(e00[i], gf, p00);
@@ -742,7 +742,7 @@ bool createRegularTwoCircleGridPeriodic (Centerline *center, GFace *gf)
   			     v2, p2,
   			     e02, pe02,
   			     e00, pe00,
-  			     e22, pe22, 
+  			     e22, pe22,
   			     quads, newv, uv, tab);
 
   printQuads(gf, uv, quads, 0);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e7ff6d65bf84ce38d6b9f283a2f8b5b3a976a123..f35fc4b03fcddf5577804280b528ef0857083fae 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -888,7 +888,7 @@ void createRegularMesh (GFace *gf,
 
  char name3[234];
   sprintf(name3,"quadParam_%d.pos", gf->tag());
-  FILE *f3 = fopen(name3,"w");
+  FILE *f3 = Fopen(name3,"w");
   fprintf(f3,"View \"%s\" {\n",name3);
 
   //create points using transfinite interpolation
@@ -1216,7 +1216,7 @@ struct  quadBlob {
     if (!CTX::instance()->mesh.saveAll) return;
     char name[234];
     sprintf(name,"blob%d_%d_%d.pos", iBlob, iter, method);
-    FILE *f = fopen(name,"w");
+    FILE *f = Fopen(name,"w");
     fprintf(f,"View \"%s\" {\n",name);
     for (unsigned int i = 0; i < quads.size(); i++){
       quads[i]->writePOS(f,true,false,false,false,false,false);
@@ -1888,7 +1888,7 @@ struct opti_data_vertex_relocation {
   }
 
   void print_cavity (char *name){
-    FILE *f = fopen(name,"w");
+    FILE *f = Fopen(name,"w");
     fprintf(f,"View \"\"{\n");
     for (unsigned int i=0;i<e.size();++i){
       MElement *el = e[i];
@@ -3586,157 +3586,157 @@ void quadsToTriangles(GFace *gf, double minqual)
 
 double Temporary::w1,Temporary::w2,Temporary::w3;
 
-std::vector<SVector3> Temporary::gradients; 	 	 
-
-Temporary::Temporary(){} 	 	 
-
-Temporary::~Temporary(){} 	 	 
-
-SVector3 Temporary::compute_gradient(MElement*element) 	 	 
-{ 	 	 
-  /*double x1,y1,z1; 	 	 
-  double x2,y2,z2; 	 	 
-  double x3,y3,z3; 	 	 
-  double x,y,z; 	 	 
-  MVertex*vertex1 = element->getVertex(0); 	 	 
-  MVertex*vertex2 = element->getVertex(1); 	 	 
-  MVertex*vertex3 = element->getVertex(2); 	 	 
-  x1 = vertex1->x(); 	 	 
-  y1 = vertex1->y(); 	 	 
-  z1 = vertex1->z(); 	 	 
-  x2 = vertex2->x(); 	 	 
-  y2 = vertex2->y(); 	 	 
-  z2 = vertex2->z(); 	 	 
-  x3 = vertex3->x(); 	 	 
-  y3 = vertex3->y(); 	 	 
-  z3 = vertex3->z(); 	 	 
-  x = (x1+x2+x3)/3.0; 	 	 
-  y = (y1+y2+y3)/3.0; 	 	 
-  z = (z1+z2+z3)/3.0;*/ 	 	 
-  return SVector3(0.0,1.0,0.0); 	 	 
-} 	 	 
-
-void Temporary::select_weights(double new_w1,double new_w2,double new_w3) 	 	 
-{ 	 	 
-  w1 = new_w1; 	 	 
-  w2 = new_w2; 	 	 
-  w3 = new_w3; 	 	 
-} 	 	 
-
-double Temporary::compute_total_cost(double f1,double f2) 	 	 
-{ 	 	 
-  double cost; 	 	 
-  cost = w1*f1 + w2*f2 + w3*f1*f2; 	 	 
-  return cost; 	 	 
-} 	 	 
-
-SVector3 Temporary::compute_normal(MElement*element) 	 	 
-{ 	 	 
-  double x1,y1,z1; 	 	 
-  double x2,y2,z2; 	 	 
-  double x3,y3,z3; 	 	 
-  double length; 	 	 
-  SVector3 vectorA; 	 	 
-  SVector3 vectorB; 	 	 
-  SVector3 normal; 	 	 
-  MVertex*vertex1 = element->getVertex(0); 	 	 
-  MVertex*vertex2 = element->getVertex(1); 	 	 
-  MVertex*vertex3 = element->getVertex(2); 	 	 
-  x1 = vertex1->x(); 	 	 
-  y1 = vertex1->y(); 	 	 
-  z1 = vertex1->z(); 	 	 
-  x2 = vertex2->x(); 	 	 
-  y2 = vertex2->y(); 	 	 
-  z2 = vertex2->z(); 	 	 
-  x3 = vertex3->x(); 	 	 
-  y3 = vertex3->y(); 	 	 
-  z3 = vertex3->z(); 	 	 
-  vectorA = SVector3(x2-x1,y2-y1,z2-z1); 	 	 
-  vectorB = SVector3(x3-x1,y3-y1,z3-z1); 	 	 
-  normal = crossprod(vectorA,vectorB); 	 	 
-  length = norm(normal); 	 	 
-  return SVector3(normal.x()/length,normal.y()/length,normal.z()/length); 	 	 
-} 	 	 
-
-SVector3 Temporary::compute_other_vector(MElement*element) 	 	 
-{ 	 	 
-  //int number; 	 	 
-  double length; 	 	 
-  SVector3 normal; 	 	 
-  SVector3 gradient; 	 	 
-  SVector3 other_vector; 	 	 
-  //number = element->getNum(); 	 	 
-  normal = Temporary::compute_normal(element); 	 	 
-  gradient = Temporary::compute_gradient(element);//gradients[number]; 	 	 
-  other_vector = crossprod(gradient,normal); 	 	 
-  length = norm(other_vector); 	 	 
-  return SVector3(other_vector.x()/length,other_vector.y()/length,other_vector.z()/length); 	 	 
-} 	 	 
-
-double Temporary::compute_alignment(const MEdge&_edge, MElement*element1, MElement*element2) 	 	 
-{ 	 	 
-  //int number; 	 	 
-  double scalar_productA,scalar_productB; 	 	 
-  double alignment; 	 	 
-  SVector3 gradient; 	 	 
-  SVector3 other_vector; 	 	 
-  SVector3 edge; 	 	 
-  MVertex*vertexA; 	 	 
-  MVertex*vertexB; 	 	 
-  //number = element1->getNum(); 	 	 
-  gradient = Temporary::compute_gradient(element1);//gradients[number]; 	 	 
-  other_vector = Temporary::compute_other_vector(element1); 	 	 
-  vertexA = _edge.getVertex(0); 	 	 
-  vertexB = _edge.getVertex(1); 	 	 
-  edge = SVector3(vertexB->x()-vertexA->x(),vertexB->y()-vertexA->y(),vertexB->z()-vertexA->z()); 	 	 
-  edge = edge * (1/norm(edge)); 	 	 
-  scalar_productA = fabs(dot(gradient,edge)); 	 	 
-  scalar_productB = fabs(dot(other_vector,edge)); 	 	 
-  alignment = std::max(scalar_productA,scalar_productB) - sqrt(2.0)/2.0; 	 	 
-  alignment = alignment/(1.0-sqrt(2.0)/2.0); 	 	 
-  return alignment; 	 	 
-} 	 	 
-
-void Temporary::read_data(std::string file_name) 	 	 
-{ 	 	 
-  #if defined(HAVE_POST) 	 	 
-  int i,j,number; 	 	 
-  double x,y,z; 	 	 
-  MElement*element; 	 	 
-  PViewData*data; 	 	 
-  PView::readMSH(file_name,-1); 	 	 
-  data = PView::list[0]->getData(); 	 	 
-  for(i = 0;i < data->getNumEntities(0);i++) 	 	 
-  { 	 	 
-    if(data->skipEntity(0,i)) continue; 	 	 
-	for(j = 0;j < data->getNumElements(0,i);j++){ 	 	 
-	  if(data->skipElement(0,i,j)) continue; 	 	 
-	  element = data->getElement(0,i,j); 	 	 
-	  number = element->getNum(); 	 	 
-	  data->getValue(0,i,j,0,x); 	 	 
-	  data->getValue(0,i,j,1,y); 	 	 
-	  data->getValue(0,i,j,2,z); 	 	 
-	  gradients[number] = SVector3(x,y,z); 	 	 
-	} 	 	 
-  } 	 	 
-  #endif 	 	 
-} 	 	 
-
-void Temporary::quadrilaterize(std::string file_name,double _w1,double _w2,double _w3) 	 	 
-{ 	 	 
-  GFace*face; 	 	 
-  GModel*model = GModel::current(); 	 	 
-  GModel::fiter iterator; 	 	 
-  gradients.resize(model->getNumMeshElements() + 1); 	 	 
-  w1 = _w1; 	 	 
-  w2 = _w2; 	 	 
-  w3 = _w3; 	 	 
-  Temporary::read_data(file_name); 	 	 
-  for(iterator = model->firstFace();iterator != model->lastFace();iterator++) 	 	 
-  { 	 	 
-    face = *iterator; 	 	 
-	_recombineIntoQuads(face,1,1); 	 	 
-  } 	 	 
+std::vector<SVector3> Temporary::gradients;
+
+Temporary::Temporary(){}
+
+Temporary::~Temporary(){}
+
+SVector3 Temporary::compute_gradient(MElement*element)
+{
+  /*double x1,y1,z1;
+  double x2,y2,z2;
+  double x3,y3,z3;
+  double x,y,z;
+  MVertex*vertex1 = element->getVertex(0);
+  MVertex*vertex2 = element->getVertex(1);
+  MVertex*vertex3 = element->getVertex(2);
+  x1 = vertex1->x();
+  y1 = vertex1->y();
+  z1 = vertex1->z();
+  x2 = vertex2->x();
+  y2 = vertex2->y();
+  z2 = vertex2->z();
+  x3 = vertex3->x();
+  y3 = vertex3->y();
+  z3 = vertex3->z();
+  x = (x1+x2+x3)/3.0;
+  y = (y1+y2+y3)/3.0;
+  z = (z1+z2+z3)/3.0;*/
+  return SVector3(0.0,1.0,0.0);
+}
+
+void Temporary::select_weights(double new_w1,double new_w2,double new_w3)
+{
+  w1 = new_w1;
+  w2 = new_w2;
+  w3 = new_w3;
+}
+
+double Temporary::compute_total_cost(double f1,double f2)
+{
+  double cost;
+  cost = w1*f1 + w2*f2 + w3*f1*f2;
+  return cost;
+}
+
+SVector3 Temporary::compute_normal(MElement*element)
+{
+  double x1,y1,z1;
+  double x2,y2,z2;
+  double x3,y3,z3;
+  double length;
+  SVector3 vectorA;
+  SVector3 vectorB;
+  SVector3 normal;
+  MVertex*vertex1 = element->getVertex(0);
+  MVertex*vertex2 = element->getVertex(1);
+  MVertex*vertex3 = element->getVertex(2);
+  x1 = vertex1->x();
+  y1 = vertex1->y();
+  z1 = vertex1->z();
+  x2 = vertex2->x();
+  y2 = vertex2->y();
+  z2 = vertex2->z();
+  x3 = vertex3->x();
+  y3 = vertex3->y();
+  z3 = vertex3->z();
+  vectorA = SVector3(x2-x1,y2-y1,z2-z1);
+  vectorB = SVector3(x3-x1,y3-y1,z3-z1);
+  normal = crossprod(vectorA,vectorB);
+  length = norm(normal);
+  return SVector3(normal.x()/length,normal.y()/length,normal.z()/length);
+}
+
+SVector3 Temporary::compute_other_vector(MElement*element)
+{
+  //int number;
+  double length;
+  SVector3 normal;
+  SVector3 gradient;
+  SVector3 other_vector;
+  //number = element->getNum();
+  normal = Temporary::compute_normal(element);
+  gradient = Temporary::compute_gradient(element);//gradients[number];
+  other_vector = crossprod(gradient,normal);
+  length = norm(other_vector);
+  return SVector3(other_vector.x()/length,other_vector.y()/length,other_vector.z()/length);
+}
+
+double Temporary::compute_alignment(const MEdge&_edge, MElement*element1, MElement*element2)
+{
+  //int number;
+  double scalar_productA,scalar_productB;
+  double alignment;
+  SVector3 gradient;
+  SVector3 other_vector;
+  SVector3 edge;
+  MVertex*vertexA;
+  MVertex*vertexB;
+  //number = element1->getNum();
+  gradient = Temporary::compute_gradient(element1);//gradients[number];
+  other_vector = Temporary::compute_other_vector(element1);
+  vertexA = _edge.getVertex(0);
+  vertexB = _edge.getVertex(1);
+  edge = SVector3(vertexB->x()-vertexA->x(),vertexB->y()-vertexA->y(),vertexB->z()-vertexA->z());
+  edge = edge * (1/norm(edge));
+  scalar_productA = fabs(dot(gradient,edge));
+  scalar_productB = fabs(dot(other_vector,edge));
+  alignment = std::max(scalar_productA,scalar_productB) - sqrt(2.0)/2.0;
+  alignment = alignment/(1.0-sqrt(2.0)/2.0);
+  return alignment;
+}
+
+void Temporary::read_data(std::string file_name)
+{
+  #if defined(HAVE_POST)
+  int i,j,number;
+  double x,y,z;
+  MElement*element;
+  PViewData*data;
+  PView::readMSH(file_name,-1);
+  data = PView::list[0]->getData();
+  for(i = 0;i < data->getNumEntities(0);i++)
+  {
+    if(data->skipEntity(0,i)) continue;
+	for(j = 0;j < data->getNumElements(0,i);j++){
+	  if(data->skipElement(0,i,j)) continue;
+	  element = data->getElement(0,i,j);
+	  number = element->getNum();
+	  data->getValue(0,i,j,0,x);
+	  data->getValue(0,i,j,1,y);
+	  data->getValue(0,i,j,2,z);
+	  gradients[number] = SVector3(x,y,z);
+	}
+  }
+  #endif
+}
+
+void Temporary::quadrilaterize(std::string file_name,double _w1,double _w2,double _w3)
+{
+  GFace*face;
+  GModel*model = GModel::current();
+  GModel::fiter iterator;
+  gradients.resize(model->getNumMeshElements() + 1);
+  w1 = _w1;
+  w2 = _w2;
+  w3 = _w3;
+  Temporary::read_data(file_name);
+  for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
+  {
+    face = *iterator;
+	_recombineIntoQuads(face,1,1);
+  }
 }
 
 /***************************************************/
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 93969c3e356496d74de6092b5da0d470f9c0def7..809307b98c31989d81d1cc93bbefd2990f50503c 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -22,6 +22,7 @@
 #include "MTetrahedron.h"
 #include "MPyramid.h"
 #include "BDS.h"
+#include "OS.h"
 #include "Context.h"
 #include "GFaceCompound.h"
 #include "meshGRegionMMG3D.h"
@@ -147,7 +148,7 @@ void printVoronoi(GRegion *gr,  std::set<SPoint3> &candidates)
   //print voronoi poles
   FILE *outfile;
   //smooth_normals *snorm = gr->model()->normals;
-  outfile = fopen("nodes.pos", "w");
+  outfile = Fopen("nodes.pos", "w");
   fprintf(outfile, "View \"Voronoi poles\" {\n");
   std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.begin();
   for(; itmap != node2Tet.end(); itmap++){
@@ -186,7 +187,7 @@ void printVoronoi(GRegion *gr,  std::set<SPoint3> &candidates)
 
   //print scalar lines
   FILE *outfile2;
-  outfile2 = fopen("edges.pos", "w");
+  outfile2 = Fopen("edges.pos", "w");
   fprintf(outfile2, "View \"Voronoi edges\" {\n");
   std::map<MFace, std::vector<MTetrahedron*> , Less_Face >::iterator itmap2 = face2Tet.begin();
   for(; itmap2 != face2Tet.end(); itmap2++){
@@ -584,7 +585,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions)
         tetrahedralize(opts, &in, &out);
         Msg::Info("%d intersecting faces have been saved into 'intersect.pos'",
                   out.numberoftrifaces);
-        FILE *fp = fopen("intersect.pos", "w");
+        FILE *fp = Fopen("intersect.pos", "w");
         if(fp){
           fprintf(fp, "View \"intersections\" {\n");
           for(int i = 0; i < out.numberoftrifaces; i++){
@@ -919,7 +920,7 @@ void meshNormalsPointOutOfTheRegion(GRegion *gr)
     }
   }
 
-  // FILE *fp = fopen("debug.pos", "w");
+  // FILE *fp = Fopen("debug.pos", "w");
   // fprintf(fp, "View \"debug\" {\n");
   // for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++)
   //   for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
diff --git a/Mesh/meshGRegionBoundaryLayers.cpp b/Mesh/meshGRegionBoundaryLayers.cpp
index 2310bcd56e94bf8a2a517ca19a30669cf3ce3173..db6f25e753d20e6bcea48ad777999610e9752247 100644
--- a/Mesh/meshGRegionBoundaryLayers.cpp
+++ b/Mesh/meshGRegionBoundaryLayers.cpp
@@ -10,7 +10,7 @@
 #include "MTetrahedron.h"
 #include "MVertex.h"
 #include "Field.h"
-
+#include "OS.h"
 
 /*
                ^  ni
@@ -189,7 +189,7 @@ BoundaryLayerColumns* buildAdditionalPoints3D (GRegion *gr)
 
   // DEBUG STUFF
 
-  FILE *f = fopen ("test3D.pos","w");
+  FILE *f = Fopen ("test3D.pos","w");
   fprintf(f,"View \"\" {\n");
   for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
     MVertex *v = *it;
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index d43bc0954d0f59032bf48f92e6648de304803ecf..7038278d3a720630044d7c2a2b6e19ffb17e228d 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -140,8 +140,8 @@ void connectTets(ITER beg, ITER end, std::set<MFace, Less_Face> *allEmbeddedFace
       for (int i = 0; i < 4; i++){
         faceXtet fxt(*beg, i);
 	// if a face is embedded, do not connect tets on both sides !
-	if (!allEmbeddedFaces || 
-	    allEmbeddedFaces->find (MFace(fxt.v[0],fxt.v[1],fxt.v[2])) == 
+	if (!allEmbeddedFaces ||
+	    allEmbeddedFaces->find (MFace(fxt.v[0],fxt.v[1],fxt.v[2])) ==
 	    allEmbeddedFaces->end()){
 	  std::set<faceXtet>::iterator found = conn.find(fxt);
 	  if (found == conn.end())
@@ -367,7 +367,7 @@ void nonrecurFindCavity(std::list<faceXtet> & shell,
 
 void printTets (const char *fn, std::list<MTet4*> &cavity, bool force = false )
 {
-  FILE *f = fopen (fn,"w");
+  FILE *f = Fopen (fn,"w");
   fprintf(f,"View \"\"{\n");
   std::list<MTet4*>::iterator ittet = cavity.begin();
   std::list<MTet4*>::iterator ittete = cavity.end();
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 77200bd1c001430177eab417af7de14c0ebd6f9c..23d41a4d3dd3b4310ee3a0a14b3155480b2fcb35 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -16,6 +16,7 @@
 #include "GFaceCompound.h"
 #include "Numeric.h"
 #include "Context.h"
+#include "OS.h"
 
 static void recur_connect(MVertex *v,
                           std::multimap<MVertex*,MEdge> &v2e,
@@ -273,7 +274,7 @@ static void printLevel(std::vector<MElement *> &elements, int recur, int region)
   }
 
   bool binary = false;
-  FILE *fp = fopen (fn, "w");
+  FILE *fp = Fopen (fn, "w");
   fprintf(fp, "$MeshFormat\n");
   fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
   fprintf(fp, "$EndMeshFormat\n");
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index a90921e3b2259e2a738f33448eea46e8060d02bc..f54e73e5d65c38b3ffd402aa587c306b81e11a7f 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -1,11 +1,20 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributor(s):
+//   Tristan Carrier Baudoin
+
 #include "GmshConfig.h"
 #include "surfaceFiller.h"
 #include "Field.h"
 #include "GModel.h"
+#include "OS.h"
 #include <queue>
 #include <stack>
 
-/// Here, we aim at producing a set of points that 
+/// Here, we aim at producing a set of points that
 /// enables to generate a nice quad mesh
 
 #if defined(HAVE_RTREE)
@@ -30,7 +39,7 @@ static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.};
 /// right spacing in the tangent plane
 
 #if defined(HAVE_RTREE)
-struct surfacePointWithExclusionRegion {  
+struct surfacePointWithExclusionRegion {
   MVertex *_v;
   SPoint2 _center;
   SPoint2 _p[4][NUMDIR];
@@ -39,7 +48,7 @@ struct surfacePointWithExclusionRegion {
   double _distanceSummed;
   /*
          + p3
-    p4   |    
+    p4   |
     +----c-----+ p2
          |
          + p1
@@ -75,7 +84,7 @@ struct surfacePointWithExclusionRegion {
     b[1] = p.y() - _q[0].y();
     sys2x2(mat, b, uv);
     //    printf("inversion 1 : %g %g \n",uv[0],uv[1]);
-    if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true; 
+    if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true;
     mat[0][0]= _q[3].x()-_q[2].x();
     mat[0][1]= _q[0].x()-_q[2].x();
     mat[1][0]= _q[3].y()-_q[2].y();
@@ -84,7 +93,7 @@ struct surfacePointWithExclusionRegion {
     b[1] = p.y() - _q[2].y();
     sys2x2(mat, b, uv);
     //    printf("inversion 2 : %g %g \n",uv[0],uv[1]);
-    if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true; 
+    if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true;
     return false;
   }
   void minmax  (double _min[2], double _max[2]) const{
@@ -100,7 +109,7 @@ struct surfacePointWithExclusionRegion {
 	    _q[1].x(),_q[1].y(),0.0,
 	    _q[2].x(),_q[2].y(),0.0,
 	    _q[3].x(),_q[3].y(),0.0,i,i,i,i);
-    
+
   }
 };
 
@@ -126,17 +135,17 @@ class compareSurfacePointWithExclusionRegionPtr
 
 bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){
   my_wrapper *w = static_cast<my_wrapper*>(point);
-  
+
   if (neighbour->inExclusionZone(w->_p)){
     w->_tooclose = true;
     return false;
   }
-  
+
   return true;
 }
 
-bool inExclusionZone (SPoint2 &p,   
-		      RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree, 
+bool inExclusionZone (SPoint2 &p,
+		      RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree,
 		      std::vector<surfacePointWithExclusionRegion*> & all ){
   // should assert that the point is inside the domain
   if (!backgroundMesh::current()->inDomain(p.x(),p.y(),0)) return true;
@@ -178,7 +187,7 @@ bool compute4neighbors (GFace *gf,   // the surface
 			MVertex *v_center, // the wertex for which we wnt to generate 4 neighbors
 			SPoint2 &midpoint,
 			bool goNonLinear, // do we compute the position in the real surface which is nonlinear
-			SPoint2 newP[4][NUMDIR], // look into other directions 
+			SPoint2 newP[4][NUMDIR], // look into other directions
 			SMetric3 &metricField, FILE *crossf = 0) // the mesh metric
 {
 
@@ -192,7 +201,7 @@ bool compute4neighbors (GFace *gf,   // the surface
 
   double L = backgroundMesh::current()->operator()(midpoint[0],midpoint[1],0.0);
   //  printf("L = %12.5E\n",L);
-  metricField = SMetric3(1./(L*L));  
+  metricField = SMetric3(1./(L*L));
   FieldManager *fields = gf->model()->getFields();
   if(fields->getBackgroundField() > 0){
     Field *f = fields->get(fields->getBackgroundField());
@@ -201,25 +210,25 @@ bool compute4neighbors (GFace *gf,   // the surface
      }
      else {
       L = (*f)(v_center->x(),v_center->y(),v_center->z(), gf);
-      metricField = SMetric3(1./(L*L));  
-     }    
+      metricField = SMetric3(1./(L*L));
+     }
   }
 
-    
+
   // get the unit normal at that point
   Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(midpoint[0],midpoint[1]));
   SVector3 s1 = der.first();
   SVector3 s2 = der.second();
   SVector3 n = crossprod(s1,s2);
   n.normalize();
-  
+
   double M = dot(s1,s1);
   double N = dot(s2,s2);
   double E = dot(s1,s2);
-  
+
 
   // compute the first fundamental form i.e. the metric tensor at the point
-  // M_{ij} = s_i \cdot s_j 
+  // M_{ij} = s_i \cdot s_j
   double metric[2][2] = {{M,E},{E,N}};
 
   //  printf("%d %g %g %g\n",gf->tag(),s1.x(),s1.y(),s1.z());
@@ -227,13 +236,13 @@ bool compute4neighbors (GFace *gf,   // the surface
   SVector3 basis_u = s1; basis_u.normalize();
   SVector3 basis_v = crossprod(n,basis_u);
 
-  for (int DIR = 0 ; DIR < NUMDIR ; DIR ++){  
+  for (int DIR = 0 ; DIR < NUMDIR ; DIR ++){
     double quadAngle  = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) + DIRS[DIR];
-    
+
     // normalize vector t1 that is tangent to gf at midpoint
     SVector3 t1 = basis_u * cos (quadAngle) + basis_v * sin (quadAngle) ;
     t1.normalize();
-    
+
     // compute the second direction t2 and normalize (t1,t2,n) is the tangent frame
     SVector3 t2 = crossprod(n,t1);
     t2.normalize();
@@ -241,7 +250,7 @@ bool compute4neighbors (GFace *gf,   // the surface
     if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),t2.x(),t2.y(),t2.z());
     if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),-t1.x(),-t1.y(),-t1.z());
     if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),-t2.x(),-t2.y(),-t2.z());
-      
+
     double size_1 = sqrt(1. / dot(t1,metricField,t1));
     double size_2 = sqrt(1. / dot(t2,metricField,t2));
 
@@ -263,14 +272,14 @@ bool compute4neighbors (GFace *gf,   // the surface
       covar2[0] = 1.0; covar2[1] = 0.0;
       singular = true;
     }
-    
+
     // transform the sizes with respect to the metric
     // consider a vector v of size 1 in the parameter plane
     // its length is sqrt (v^T M v) --> if I want a real size
     // of size1 in direction v, it should be sqrt(v^T M v) * size1
     double l1 = sqrt(covar1[0]*covar1[0]+covar1[1]*covar1[1]);
     double l2 = sqrt(covar2[0]*covar2[0]+covar2[1]*covar2[1]);
-  
+
     covar1[0] /= l1;covar1[1] /= l1;
     covar2[0] /= l2;covar2[1] /= l2;
 
@@ -278,7 +287,7 @@ bool compute4neighbors (GFace *gf,   // the surface
 						  2*E*covar1[1]*covar1[0]+
 						  N*covar1[1]*covar1[1]);
     double size_param_2  = size_2 / sqrt (  M*covar2[0]*covar2[0]+
-						  2*E*covar2[1]*covar2[0]+					    
+						  2*E*covar2[1]*covar2[0]+
 						  N*covar2[1]*covar2[1]);
     if (singular){
       size_param_1 = size_param_2 = std::min (size_param_1,size_param_2);
@@ -286,7 +295,7 @@ bool compute4neighbors (GFace *gf,   // the surface
 
     //    printf("%12.5E %12.5E\n", size_param_1, size_param_2);
 
-    
+
     //    if (v_center->onWhat() != gf && gf->tag() == 3)
     //      printf("M = (%g %g %g) L = %g %g LP = %g %g\n",metricField(0,0),metricField(1,1),metricField(0,1),l1,l2,size_param_1,size_param_2);
     //if (l1 == 0.0 || l2 == 0.0) printf("bouuuuuuuuuuuuh %g %g %g %g --- %g %g %g %g %g %g\n",l1,l2,t1.norm(),t2.norm(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z());
@@ -294,7 +303,7 @@ bool compute4neighbors (GFace *gf,   // the surface
     /*    printf("%12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %g %g %g %g %g %g %g %g %g %g %g\n",
 	   M*covar1[0]*covar1[0]+
 	   2*E*covar1[1]*covar1[0]+
-	   N*covar1[1]*covar1[1], 
+	   N*covar1[1]*covar1[1],
 	   M*covar2[0]*covar2[0]+
 	   2*E*covar2[1]*covar2[0]+
 	   N*covar2[1]*covar2[1]
@@ -323,13 +332,13 @@ bool compute4neighbors (GFace *gf,   // the surface
     // surface
     double ERR[4];
     for (int i=0;i<4;i++){                                              //
-      //      if (newPoint[i][0] < rangeU.low())newPoint[i][0] = rangeU.low(); 
-      //      if (newPoint[i][0] > rangeU.high())newPoint[i][0] = rangeU.high(); 
-      //      if (newPoint[i][1] < rangeV.low())newPoint[i][1] = rangeV.low(); 
-      //      if (newPoint[i][1] > rangeV.high())newPoint[i][1] = rangeV.high(); 
+      //      if (newPoint[i][0] < rangeU.low())newPoint[i][0] = rangeU.low();
+      //      if (newPoint[i][0] > rangeU.high())newPoint[i][0] = rangeU.high();
+      //      if (newPoint[i][1] < rangeV.low())newPoint[i][1] = rangeV.low();
+      //      if (newPoint[i][1] > rangeV.high())newPoint[i][1] = rangeV.high();
       GPoint pp = gf->point(SPoint2(newPoint[i][0], newPoint[i][1]));
-      double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) + 
-		       (pp.y() - v_center->y())*(pp.y() - v_center->y()) + 
+      double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
+		       (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
 		       (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
       ERR[i] = 100*fabs(D-L)/(D+L);
       //      printf("L = %12.5E D = %12.5E ERR = %12.5E\n",L,D,100*fabs(D-L)/(D+L));
@@ -337,31 +346,31 @@ bool compute4neighbors (GFace *gf,   // the surface
 
     if (1 && goNonLinear){//---------------------------------------------------//
       surfaceFunctorGFace ss (gf);                                        //
-      SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1*(1.0),t2*(1.0)};                     //      
+      SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1*(1.0),t2*(1.0)};                     //
       for (int i=0;i<4;i++){                                              //
 	if (ERR[i] > 12){
 	  double uvt[3] = {newPoint[i][0],newPoint[i][1],0.0};              //
 	  //	  printf("Intersecting with circle N = %g %g %g dir = %g %g %g R = %g p = %g %g %g\n",n.x(),n.y(),n.z(),dirs[i].x(),dirs[i].y(),dirs[i].z(),L,v_center->x(),v_center->y(),v_center->z());
-	  curveFunctorCircle cf (dirs[i],n,				 
+	  curveFunctorCircle cf (dirs[i],n,
 				 SVector3(v_center->x(),v_center->y(),v_center->z()),
 				 L);
 	  if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-3)){          //
 	    GPoint pp = gf->point(SPoint2(uvt[0],uvt[1]));
-	    double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) + 
-			     (pp.y() - v_center->y())*(pp.y() - v_center->y()) + 
+	    double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
+			     (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
 			     (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
-	    double DP = sqrt ((newPoint[i][0]-uvt[0])*(newPoint[i][0]-uvt[0]) + 
+	    double DP = sqrt ((newPoint[i][0]-uvt[0])*(newPoint[i][0]-uvt[0]) +
 			      (newPoint[i][1]-uvt[1])*(newPoint[i][1]-uvt[1]));
 	    double newErr = 100*fabs(D-L)/(D+L);
 	    //	    if (v_center->onWhat() != gf && gf->tag() == 3){
 	    //	      crossField2d::normalizeAngle (uvt[2]);
-	    //	      printf("INTERSECT angle = %g DP %g\n",uvt[2],DP); 
+	    //	      printf("INTERSECT angle = %g DP %g\n",uvt[2],DP);
 	    //	    }
 	    if (newErr < 1 && DP < .1){
 	      //	      printf("%12.5E vs %12.5E : %12.5E  %12.5E vs %12.5E  %12.5E \n",ERR[i],newErr,newPoint[i][0],newPoint[i][1],uvt[0],uvt[1]);
 	      newPoint[i][0] = uvt[0];                                        //
 	      newPoint[i][1] = uvt[1];                                        //
-	    }                                                                 //	    
+	    }                                                                 //
 	    //	    printf("OK\n");
 	  }
 	  else{
@@ -371,14 +380,14 @@ bool compute4neighbors (GFace *gf,   // the surface
 	}
       }                                                                   //
     } /// end non linear -------------------------------------------------//
-    
+
     // return the four new vertices
     for (int i=0;i<4;i++){
       newP[i][DIR] = SPoint2(newPoint[i][0],newPoint[i][1]);
     }
   }
   return true;
-} 
+}
 #endif
 
 // fills a surface with points in order to build a nice
@@ -388,8 +397,8 @@ void packingOfParallelograms(GFace* gf,  std::vector<MVertex*> &packed, std::vec
 
   const bool goNonLinear = true;
 
-  //  FILE *f = fopen ("parallelograms.pos","w"); 
-  
+  //  FILE *f = Fopen ("parallelograms.pos","w");
+
   // get all the boundary vertices
   std::set<MVertex*> bnd_vertices;
   for(unsigned int i=0;i<gf->getNumMeshElements();i++){
@@ -400,7 +409,7 @@ void packingOfParallelograms(GFace* gf,  std::vector<MVertex*> &packed, std::vec
     }
   }
 
-  // put boundary vertices in a fifo queue  
+  // put boundary vertices in a fifo queue
   // std::queue<surfacePointWithExclusionRegion*> fifo;
   std::set<surfacePointWithExclusionRegion*,  compareSurfacePointWithExclusionRegionPtr> fifo;
   std::vector<surfacePointWithExclusionRegion*> vertices;
@@ -411,16 +420,16 @@ void packingOfParallelograms(GFace* gf,  std::vector<MVertex*> &packed, std::vec
   std::set<MVertex*>::iterator it =  bnd_vertices.begin() ;
 
   char NAME[345]; sprintf(NAME,"crossReal%d.pos",gf->tag());
-  FILE *crossf = fopen (NAME,"w");
+  FILE *crossf = Fopen (NAME,"w");
   if (crossf)fprintf(crossf,"View \"\"{\n");
   for (; it !=  bnd_vertices.end() ; ++it){
     SPoint2 midpoint;
     compute4neighbors (gf, *it, midpoint, goNonLinear, newp, metricField,crossf);
-    surfacePointWithExclusionRegion *sp = 
-      new surfacePointWithExclusionRegion (*it, newp, midpoint,metricField);    
-    //    fifo.push(sp); 
-    fifo.insert(sp); 
-    vertices.push_back(sp); 
+    surfacePointWithExclusionRegion *sp =
+      new surfacePointWithExclusionRegion (*it, newp, midpoint,metricField);
+    //    fifo.push(sp);
+    fifo.insert(sp);
+    vertices.push_back(sp);
     double _min[2],_max[2];
     sp->minmax(_min,_max);
     //    printf("%g %g .. %g %g\n",_min[0],_min[1],_max[0],_max[1]);
@@ -442,29 +451,29 @@ void packingOfParallelograms(GFace* gf,  std::vector<MVertex*> &packed, std::vec
       int countOK = 0;
       for (int i=0;i<4;i++){
 	//	printf("i = %d %12.5E %12.5E \n",i,parent._p[i][dir].x(),parent._p[i][dir].y());
-		
-	//	if (!w._tooclose){	
-	if (!inExclusionZone (parent->_p[i][dir], rtree, vertices) ){	
+
+	//	if (!w._tooclose){
+	if (!inExclusionZone (parent->_p[i][dir], rtree, vertices) ){
 	  countOK++;
 	  GPoint gp = gf->point(parent->_p[i][dir]);
 	  MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
 	  //	  	printf(" %g %g %g %g\n",parent._center.x(),parent._center.y(),gp.u(),gp.v());
 	  SPoint2 midpoint;
 	  compute4neighbors (gf, v, midpoint, goNonLinear, newp, metricField,crossf);
-	  surfacePointWithExclusionRegion *sp = 
-	    new surfacePointWithExclusionRegion (v, newp, midpoint, metricField, parent);    
-	  //	  fifo.push(sp); 
-	  fifo.insert(sp); 
-	  vertices.push_back(sp); 
+	  surfacePointWithExclusionRegion *sp =
+	    new surfacePointWithExclusionRegion (v, newp, midpoint, metricField, parent);
+	  //	  fifo.push(sp);
+	  fifo.insert(sp);
+	  vertices.push_back(sp);
 	  double _min[2],_max[2];
 	  sp->minmax(_min,_max);
-	  rtree.Insert(_min,_max,sp);    
+	  rtree.Insert(_min,_max,sp);
 	}
       }
       if (countOK)break;
     }
     //    printf("%d\n",vertices.size());
-  }  
+  }
   if (crossf){
     fprintf(crossf,"};\n");
     fclose (crossf);
@@ -474,10 +483,10 @@ void packingOfParallelograms(GFace* gf,  std::vector<MVertex*> &packed, std::vec
   // add the vertices as additional vertices in the
   // surface mesh
   char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
-  FILE *f = fopen(ccc,"w");
+  FILE *f = Fopen(ccc,"w");
   fprintf(f,"View \"\"{\n");
   for (unsigned int i=0;i<vertices.size();i++){
-    //    if(vertices[i]->_v->onWhat() != gf) 
+    //    if(vertices[i]->_v->onWhat() != gf)
       vertices[i]->print(f,i);
     if(vertices[i]->_v->onWhat() == gf) {
       packed.push_back(vertices[i]->_v);
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 46463cf704729206b7a1474c5d9b10a73494dc96..57c4cec64c5b2ec8eae79a6918be79f3f85ab472 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -4455,7 +4455,7 @@ yyreduce:
 #line 220 "Gmsh.y"
     {
       std::string tmp = FixRelativePath(gmsh_yyname, (yyvsp[(6) - (7)].c));
-      FILE *fp = fopen(tmp.c_str(), (yyvsp[(5) - (7)].c));
+      FILE *fp = Fopen(tmp.c_str(), (yyvsp[(5) - (7)].c));
       if(!fp){
 	yymsg(0, "Unable to open file '%s'", tmp.c_str());
       }
@@ -4511,7 +4511,7 @@ yyreduce:
 	yymsg(0, "%d extra argument%s in Printf", i, (i > 1) ? "s" : "");
       else{
         std::string tmp = FixRelativePath(gmsh_yyname, (yyvsp[(8) - (9)].c));
-	FILE *fp = fopen(tmp.c_str(), (yyvsp[(7) - (9)].c));
+	FILE *fp = Fopen(tmp.c_str(), (yyvsp[(7) - (9)].c));
 	if(!fp){
 	  yymsg(0, "Unable to open file '%s'", tmp.c_str());
 	}
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 081cf9ed1ca66f564f5883cd5310c63faf152b57..cced76e308f0feb2356e57027b3376ded766cdb7 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -219,7 +219,7 @@ Printf :
   | tPrintf '(' tBIGSTR ')' SendToFile StringExprVar tEND
     {
       std::string tmp = FixRelativePath(gmsh_yyname, $6);
-      FILE *fp = fopen(tmp.c_str(), $5);
+      FILE *fp = Fopen(tmp.c_str(), $5);
       if(!fp){
 	yymsg(0, "Unable to open file '%s'", tmp.c_str());
       }
@@ -266,7 +266,7 @@ Printf :
 	yymsg(0, "%d extra argument%s in Printf", i, (i > 1) ? "s" : "");
       else{
         std::string tmp = FixRelativePath(gmsh_yyname, $8);
-	FILE *fp = fopen(tmp.c_str(), $7);
+	FILE *fp = Fopen(tmp.c_str(), $7);
 	if(!fp){
 	  yymsg(0, "Unable to open file '%s'", tmp.c_str());
 	}
diff --git a/Plugin/Bubbles.cpp b/Plugin/Bubbles.cpp
index 419f9313f46dbc7b3445c31861c9ee5826357e1c..2155e85f992f00af79f4d4e580bae0602936aea2 100644
--- a/Plugin/Bubbles.cpp
+++ b/Plugin/Bubbles.cpp
@@ -9,6 +9,7 @@
 #include "MTriangle.h"
 #include "Numeric.h"
 #include "Bubbles.h"
+#include "OS.h"
 
 StringXNumber BubblesOptions_Number[] = {
   {GMSH_FULLRC, "ShrinkFactor", NULL, 0.},
@@ -60,8 +61,8 @@ static double myangle(double c[3], double p[3])
 {
   double v[3] = {1, 0, 0};
   double n[3] = {0, 0, 1};
-  if(fabs(c[0] - p[0]) < 1e-12 && 
-     fabs(c[1] - p[1]) < 1e-12 && 
+  if(fabs(c[0] - p[0]) < 1e-12 &&
+     fabs(c[1] - p[1]) < 1e-12 &&
      fabs(c[2] - p[2]) < 1e-12)
     return 0.;
   return angle_plan(c, v, p, n);
@@ -88,12 +89,12 @@ PView *GMSH_BubblesPlugin::execute(PView *v)
   double shrink = (double)BubblesOptions_Number[0].def;
   std::string fileName = BubblesOptions_String[0].def;
 
-  FILE *fp = fopen(fileName.c_str(), "w");
+  FILE *fp = Fopen(fileName.c_str(), "w");
   if(!fp){
     Msg::Error("Could not open output file '%s'", fileName.c_str());
     return v;
   }
-  
+
   GModel *m = GModel::current();
 
   int p = m->getMaxElementaryNumber(0) + 1;
@@ -131,7 +132,7 @@ PView *GMSH_BubblesPlugin::execute(PView *v)
       if(v->onWhat() && v->onWhat()->dim() < 2)
         it->second.push_back(SPoint3(it->first->x(), it->first->y(), it->first->z()));
     }
- 
+
     for(std::map<MVertex*, std::vector<SPoint3> >::iterator it = v2t.begin();
         it != v2t.end(); it++){
       if(it->second.size() > 2){
@@ -145,8 +146,8 @@ PView *GMSH_BubblesPlugin::execute(PView *v)
         // shrink cells
         if(shrink){
           for(unsigned int i = 0; i < it->second.size(); i++){
-            double dir[3] = {it->second[i].x() - bc.x(), 
-                             it->second[i].y() - bc.y(), 
+            double dir[3] = {it->second[i].x() - bc.x(),
+                             it->second[i].y() - bc.y(),
                              it->second[i].z() - bc.z()};
             it->second[i][0] -= shrink * dir[0];
             it->second[i][1] -= shrink * dir[1];
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 9eea7c6a96fede61b339fa2111db9e287bdcf8f0..a652bfe20b769bf537ad72c91f3b10d833c4dd78 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -7,6 +7,7 @@
 #include "Gmsh.h"
 #include "GmshConfig.h"
 #include "GModel.h"
+#include "OS.h"
 #include "Distance.h"
 #include "distanceTerm.h"
 #include "Context.h"
@@ -108,7 +109,7 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
   }
 
   Msg::Info("Writing %s", _fileName.c_str());
-  FILE *fName = fopen(_fileName.c_str(),"w");
+  FILE *fName = Fopen(_fileName.c_str(),"w");
   fprintf(fName, "View \"distance \"{\n");
 
   for (unsigned int ii=0; ii<_entities.size(); ii++) {
@@ -500,7 +501,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
     data2->setName("ortogonal field");
 
     Msg::Info("Writing  orthogonal.pos");
-    FILE * f5 = fopen("orthogonal.pos","w");
+    FILE * f5 = Fopen("orthogonal.pos","w");
     fprintf(f5,"View \"orthogonal\"{\n");
     for (std::vector<MElement* >::iterator it = allElems.begin();
          it != allElems.end(); it++){
diff --git a/Plugin/NearToFarField.cpp b/Plugin/NearToFarField.cpp
index 771728912c53582a0ce53e5463809ad3d32e68cb..b352f7b9df5a38705449685aa592cf8fb2013f9f 100644
--- a/Plugin/NearToFarField.cpp
+++ b/Plugin/NearToFarField.cpp
@@ -9,6 +9,7 @@
 
 #include <complex>
 #include "NearToFarField.h"
+#include "OS.h"
 
 StringXNumber NearToFarFieldOptions_Number[] = {
   {GMSH_FULLRC, "Wavenumber",       NULL, 1.},
@@ -460,7 +461,7 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
   }
 
   if(_outFile.size()){
-    FILE *fp = fopen(_outFile.c_str(), "w");
+    FILE *fp = Fopen(_outFile.c_str(), "w");
     if(fp){
       printVector(fp, "phi", phi);
       printVector(fp, "theta", theta);
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index d177ce0c7b6dde99684a0a9415a5c43870855925..6d79ff9ad3394ad3b7f9cbe0946e0613d4bfe03c 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -10,6 +10,7 @@
 #include "MElement.h"
 #include "Numeric.h"
 #include "StringUtils.h"
+#include "OS.h"
 
 bool PViewDataGModel::addData(GModel *model, std::map<int, std::vector<double> > &data,
                               int step, double time, int partition, int numComp)
@@ -149,7 +150,7 @@ bool PViewDataGModel::writeMSH(const std::string &fileName, double version, bool
     if(!model->writeMSH(fileName, version, binary, false, false, 1.0, 0,
                         0, multipleView)) return false;
     // append data
-    fp = fopen(fileName.c_str(), binary ? "ab" : "a");
+    fp = Fopen(fileName.c_str(), binary ? "ab" : "a");
     if(!fp){
       Msg::Error("Unable to open file '%s'", fileName.c_str());
       return false;
@@ -157,14 +158,14 @@ bool PViewDataGModel::writeMSH(const std::string &fileName, double version, bool
   }
   else{
     if(multipleView){
-      fp = fopen(fileName.c_str(), binary ? "ab" : "a");
+      fp = Fopen(fileName.c_str(), binary ? "ab" : "a");
       if(!fp){
         Msg::Error("Unable to open file '%s'", fileName.c_str());
         return false;
       }
     }
     else{
-      fp = fopen(fileName.c_str(), binary ? "wb" : "w");
+      fp = Fopen(fileName.c_str(), binary ? "wb" : "w");
       if(!fp){
         Msg::Error("Unable to open file '%s'", fileName.c_str());
         return false;
diff --git a/Post/PViewDataIO.cpp b/Post/PViewDataIO.cpp
index 02d9da6a7e6a07a8e93ca58047b204dd14d86f11..6b571abff553716c01f5dabca87ca0c6763c11b6 100644
--- a/Post/PViewDataIO.cpp
+++ b/Post/PViewDataIO.cpp
@@ -10,10 +10,11 @@
 #include "Numeric.h"
 #include "PViewData.h"
 #include "adaptiveData.h"
+#include "OS.h"
 
 bool PViewData::writeSTL(const std::string &fileName)
 {
-  FILE *fp = fopen(fileName.c_str(), "w");
+  FILE *fp = Fopen(fileName.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", fileName.c_str());
     return false;
@@ -73,7 +74,7 @@ bool PViewData::writeSTL(const std::string &fileName)
 
 bool PViewData::writeTXT(const std::string &fileName)
 {
-  FILE *fp = fopen(fileName.c_str(), "w");
+  FILE *fp = Fopen(fileName.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", fileName.c_str());
     return false;
@@ -119,7 +120,7 @@ bool PViewData::writePOS(const std::string &fileName, bool binary, bool parsed,
   if(binary || !parsed)
     Msg::Warning("Only parsed .pos files can be exported for this view type");
 
-  FILE *fp = fopen(fileName.c_str(), append ? "a" : "w");
+  FILE *fp = Fopen(fileName.c_str(), append ? "a" : "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", fileName.c_str());
     return false;
diff --git a/Post/PViewDataListIO.cpp b/Post/PViewDataListIO.cpp
index 442c23e2a33fed66ddf8285e34d32b24f459fc0a..8ba69355a721d32e9003a32accc564374826de56 100644
--- a/Post/PViewDataListIO.cpp
+++ b/Post/PViewDataListIO.cpp
@@ -14,6 +14,7 @@
 #include "MVertexPositionSet.h"
 #include "Context.h"
 #include "adaptiveData.h"
+#include "OS.h"
 
 static void dVecRead(std::vector<double> &v, int n, FILE *fp,
                      bool binary, int swap)
@@ -356,7 +357,7 @@ bool PViewDataList::writePOS(const std::string &fileName, bool binary, bool pars
     return false;
   }
 
-  FILE *fp = fopen(fileName.c_str(),
+  FILE *fp = Fopen(fileName.c_str(),
                    append ? (binary ? "ab" : "a") : (binary ? "wb" : "w"));
   if(!fp){
     Msg::Error("Unable to open file '%s'", fileName.c_str());
@@ -533,7 +534,7 @@ bool PViewDataList::writeMSH(const std::string &fileName, double version, bool b
     return _adaptive->getData()->writeMSH(fileName, version, binary);
   }
 
-  FILE *fp = fopen(fileName.c_str(), "w");
+  FILE *fp = Fopen(fileName.c_str(), "w");
   if(!fp){
     Msg::Error("Unable to open file '%s'", fileName.c_str());
     return false;
@@ -661,7 +662,7 @@ void PViewDataList::getListPointers(int N[24], std::vector<double> *V[24])
     std::vector<double> *list = 0;
     int *nbe = 0, nbc, nbn;
     _getRawData(i, &list, &nbe, &nbc, &nbn);
-    N[i] = *nbe; 
+    N[i] = *nbe;
     V[i] = list; // copy pointer only
   }
 }
diff --git a/Post/PViewIO.cpp b/Post/PViewIO.cpp
index 36520e3d951cc0530f94a729b16a3d05834ded8c..df3e0edaeeb3b0a63efa516652a5253fa9959e28 100644
--- a/Post/PViewIO.cpp
+++ b/Post/PViewIO.cpp
@@ -10,10 +10,11 @@
 #include "PViewDataGModel.h"
 #include "StringUtils.h"
 #include "Context.h"
+#include "OS.h"
 
 bool PView::readPOS(const std::string &fileName, int fileIndex)
 {
-  FILE *fp = fopen(fileName.c_str(), "rb");
+  FILE *fp = Fopen(fileName.c_str(), "rb");
   if(!fp){
     Msg::Error("Unable to open file '%s'", fileName.c_str());
     return false;
@@ -83,7 +84,7 @@ bool PView::readPOS(const std::string &fileName, int fileIndex)
 
 bool PView::readMSH(const std::string &fileName, int fileIndex)
 {
-  FILE *fp = fopen(fileName.c_str(), "rb");
+  FILE *fp = Fopen(fileName.c_str(), "rb");
   if(!fp){
     Msg::Error("Unable to open file '%s'", fileName.c_str());
     return false;
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 382994776ea55d09ab8c4e7621b6b9cfeffef458..6b8d4369f622320ae28079759c027531842b5d76 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -12,6 +12,7 @@
 #include "linearSystemFull.h"
 #include "Numeric.h"
 #include "GModel.h"
+#include "OS.h"
 #include "functionSpace.h"
 #include "terms.h"
 #include "solverAlgorithms.h"
@@ -70,8 +71,7 @@ void elasticitySolver::setMesh(const std::string &meshFileName)
 
 void elasticitySolver::exportKb()
 {
-
-  FILE *f = fopen ( "K.txt", "w" );
+  FILE *f = Fopen ( "K.txt", "w" );
   double valeur;
   std::string sysname = "A";
   for ( int i = 0 ; i < pAssembler->sizeOfR() ; i ++ )
@@ -86,7 +86,7 @@ void elasticitySolver::exportKb()
 
   fclose ( f );
 
-  f = fopen ( "b.txt", "w" );
+  f = Fopen ( "b.txt", "w" );
   for ( int i = 0 ; i < pAssembler->sizeOfR() ; i ++ )
   {
     pAssembler->getLinearSystem ( sysname )->getFromRightHandSide ( i,valeur );
@@ -94,7 +94,6 @@ void elasticitySolver::exportKb()
   }
 
   fclose ( f );
-
 }
 
 void elasticitySolver::solve()
@@ -139,14 +138,15 @@ void elasticitySolver::postSolve()
     SolverField<SVector3> Field(pAssembler, LagSpace);
     IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
     BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
-    Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
+    Assemble(Elastic_Energy_Term, elasticFields[i].g->begin(),
+             elasticFields[i].g->end(), Integ_Bulk, energ);
   }
   printf("elastic energy=%f\n",energ);
 }
 
 void elasticitySolver::readInputFile(const std::string &fn)
 {
-  FILE *f = fopen(fn.c_str(), "r");
+  FILE *f = Fopen(fn.c_str(), "r");
   char what[256];
   while(!feof(f)){
     if(fscanf(f, "%s", what) != 1) return;
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 18fd5496748a8b577a2ba3101339778eeacbc38b..b7f0dcb4fd76190b7ba2a566390a2fc4391c947b 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -631,7 +631,7 @@ static void printCut(std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
 {
    printf("Writing points.pos \n");
    std::map<MEdge,MVertex*,Less_Edge>::iterator ite = cutEdges.begin();
-   FILE *f1 = fopen("points.pos","w");
+   FILE *f1 = Fopen("points.pos","w");
    fprintf(f1,"View\"\"{\n");
    for ( ; ite != cutEdges.end();++ite){
      fprintf(f1,"SP(%g,%g,%g){1.0};\n",ite->second->x(),ite->second->y(),ite->second->z());
@@ -645,7 +645,7 @@ static void printCut(std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
 
    printf("Writing edges.pos \n");
    std::set<MEdge,Less_Edge>::iterator itc = theCut.begin();
-   FILE *f2 = fopen("edges.pos","w");
+   FILE *f2 = Fopen("edges.pos","w");
    fprintf(f2,"View\"\"{\n");
    for ( ; itc != theCut.end();++itc){
      fprintf(f2,"SL(%g,%g,%g,%g,%g,%g){1.0,1.0};\n",itc->getVertex(0)->x(),
@@ -670,7 +670,7 @@ static void printLevel(const char* fn,
       vs.insert(elements[i]->getVertex(j));
 
   bool binary = false;
-  FILE *fp = fopen (fn, "w");
+  FILE *fp = Fopen (fn, "w");
   fprintf(fp, "$MeshFormat\n");
   fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
   fprintf(fp, "$EndMeshFormat\n");