diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 0aee27331f20821b99417905302e90d45cd52bb0..f9c3a9c700e9c297c76476f969d3b496510c565d 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -774,12 +774,7 @@ StringXNumber GeneralOptions_Number[] = {
   { F|S, "SystemMenuBar" , opt_general_system_menu_bar , 1. , 
     "Use the system menu bar on Mac OS X?" }, 
 
-  { F|O, "Terminal" , opt_general_terminal , 
-#if defined(HAVE_FLTK)
-    0. ,
-#else                 
-    1. ,
-#endif
+  { F|O, "Terminal" , opt_general_terminal , 0. ,
     "Should information be printed on the terminal (if available)?" },
   { F|O, "Tooltips" , opt_general_tooltips , 1. ,
     "Show tooltips in the user interface" },
diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index c6eda17d78b925daa299099ee961d56956f26fb0..a13fffa82a29233ff13a450a0e87ddbcbe86968e 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -1,4 +1,4 @@
-// $Id: Gmsh.cpp,v 1.6 2008-05-04 08:31:11 geuzaine Exp $
+// $Id: Gmsh.cpp,v 1.7 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -68,6 +68,12 @@ int GmshInitialize(int argc, char **argv)
   return 1;
 }
 
+int GmshSetMessageHandler(GmshMessage *callback)
+{
+  Msg::SetCallback(callback);
+  return 1;
+}
+
 int GmshFinalize()
 {
   return 1;
diff --git a/Common/Gmsh.h b/Common/Gmsh.h
index 2a9178c7c597a980f90816d870b560380a175c2d..fe50f022e43dfa5e552c9581116a9947eb350842 100644
--- a/Common/Gmsh.h
+++ b/Common/Gmsh.h
@@ -19,7 +19,17 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
+#include <string>
+
+class GmshMessage{
+ public:
+  GmshMessage(){}
+  virtual ~GmshMessage(){}
+  virtual void operator()(std::string level, std::string message){}
+};
+
 int GmshInitialize(int argc=0, char **argv=0);
+int GmshSetMessageHandler(GmshMessage *callback);
 int GmshFinalize();
 int GmshBatch();
 
diff --git a/Common/GmshMatrix.h b/Common/GmshMatrix.h
index 18fcff0e7cb321d3ea8a2429a7a60478a1df8bfa..f460b383ed9568c5a13ddf41a24f031961cff3b9 100644
--- a/Common/GmshMatrix.h
+++ b/Common/GmshMatrix.h
@@ -19,14 +19,12 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
-#include <assert.h>
-
 template <class SCALAR>
 class Gmsh_Vector
 {
-private:
+ private:
   int r;
-public:
+ public:
   inline int size() const { return r; }
   SCALAR *data;
   ~Gmsh_Vector() { delete [] data; }
@@ -63,9 +61,9 @@ public:
 template <class SCALAR>
 class Gmsh_Matrix
 {
-private:
+ private:
   int r, c;
-public:
+ public:
   inline int size1() const { return r; }
   inline int size2() const { return c; }
   SCALAR *data;
@@ -95,39 +93,46 @@ public:
   }
   inline void mult(const Gmsh_Matrix<SCALAR> &x, const Gmsh_Matrix<SCALAR> &b)
   {
-    throw;
-  }
-  inline void blas_dgemm(const Gmsh_Matrix<SCALAR>& x, const Gmsh_Matrix<SCALAR>& b, 
-			 const double c_a = 1.0, const double c_b = 1.0){
-    throw;
+    for(int i = 0; i < b.size1(); i++)
+      for(int j = 0; j < b.size2(); j++)
+	for(int k = 0; k < size2(); k++)
+	  b(i, j) += (*this)(i, k) * x(k, j);
   }
   inline void mult(const Gmsh_Vector<SCALAR> &x, Gmsh_Vector<SCALAR> &b)
   {
-    throw;
+    for(int i = 0; i < b.size(); i++)
+      for(int j = 0; j < b.size(); j++)
+	b(i) += (*this)(i, j) * x(j);
   }
   inline void set_all(const double &m) 
   {
-    throw;
+    for(int i = 0; i < r * c; i++) data[i] = m;
   }
   inline void lu_solve(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
   {
-    throw;
+    // FIXME: not implemented
+    result.scale(0);
   }
   Gmsh_Matrix cofactor(int i, int j) const 
   {
-    throw;
+    // FIXME: not implemented
+    Gmsh_Matrix cof(size1() - 1, size2() - 1);
+    return cof;
   }
   inline void invert()
   {
-    throw;
+    // FIXME: not implemented
   }
   double determinant() const 
   {
-    throw;
+    // FIXME: not implemented
+    return 0.;
   }
   inline Gmsh_Matrix touchSubmatrix(int i0, int ni, int j0, int nj) 
   {
-    throw;
+    // FIXME: not implemented
+    Gmsh_Matrix subm(ni, nj);
+    return subm;
   }  
   inline void scale(const SCALAR s)
   {
@@ -135,15 +140,17 @@ public:
   }
   inline void add(const double &a) 
   {
-    throw;
+    for (int i = 0; i < r * c; ++i) data[i] += a;
   }
   inline void add(const Gmsh_Matrix &m) 
   {
-    throw;
+    for(int i = 0; i < size1(); i++)
+      for(int j = 0; j < size2(); j++)
+	(*this)(i, j) += m(i, j);
   }
 };
 
-#ifdef HAVE_GSL
+#if defined(HAVE_GSL)
 
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
@@ -152,9 +159,9 @@ public:
 
 class GSL_Vector
 {
-private:
+ private:
   int r;
-public:
+ public:
   inline int size() const { return r; }
   gsl_vector *data;
   ~GSL_Vector() { gsl_vector_free(data); }
@@ -176,7 +183,7 @@ public:
     return *gsl_vector_ptr(data, i);
   }
   inline double norm()
-{
+  {
     return gsl_blas_dnrm2(data);
   }
   inline void scale(const double &y)
@@ -188,9 +195,9 @@ public:
 
 class GSL_Matrix
 {
-private:
+ private:
   gsl_matrix_view view;
-public:
+ public:
   inline int size1() const { return data->size1; }
   inline int size2() const { return data->size2; }
   gsl_matrix *data;
@@ -229,11 +236,6 @@ public:
   {
     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
   }
-  inline void blas_dgemm(const GSL_Matrix & x, const GSL_Matrix& b, 
-			 const double c_a = 1.0, const double c_b = 1.0)
-  {      
-    gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, c_a, x.data, b.data, c_b, data);
-  }
   inline void set_all(const double &m) 
   {
     gsl_matrix_set_all(data, m);
@@ -323,7 +325,7 @@ public:
   }
   inline void add(const GSL_Matrix &m) 
   {
-    gsl_matrix_add (data, m.data);
+    gsl_matrix_add(data, m.data);
   }
 };
 
diff --git a/Common/Main.cpp b/Common/Main.cpp
index 24f2b4a54f73056dd38291b512701ef657cc805d..eefd628f25801423860b039d24c8fb7339fec5eb 100644
--- a/Common/Main.cpp
+++ b/Common/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.3 2008-07-02 17:40:56 geuzaine Exp $
+// $Id: Main.cpp,v 1.4 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -41,6 +41,9 @@ int main(int argc, char *argv[])
   }
 
   GmshInitialize(argc, argv);
+  // force this even if the options say it ain't so
+  CTX.terminal = 1; 
+
   new GModel;
   GmshBatch();
   GmshFinalize();
diff --git a/Common/Makefile b/Common/Makefile
index 015f1f3b22d92169c7fb5e86b168f2463a340572..8b56ae3892218a2cd006c3a9d5420f2a5a7a56b3 100644
--- a/Common/Makefile
+++ b/Common/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.179 2008-06-07 17:22:37 geuzaine Exp $
+# $Id: Makefile,v 1.180 2008-07-03 17:06:01 geuzaine Exp $
 #
 # Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 #
@@ -126,9 +126,9 @@ Gmsh.o: Gmsh.cpp GmshDefines.h ../Geo/GModel.h ../Geo/GVertex.h \
   ../Common/Message.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Common/GmshMatrix.h
 OS.o: OS.cpp Message.h
-Message.o: Message.cpp Message.h Options.h ../Post/ColorTable.h Context.h \
-  OS.h ../Fltk/GUI.h ../Fltk/Opengl_Window.h ../Fltk/Colorbar_Window.h \
-  ../Common/GmshUI.h ../Fltk/Popup_Button.h \
+Message.o: Message.cpp Message.h Gmsh.h Options.h ../Post/ColorTable.h \
+  Context.h OS.h ../Fltk/GUI.h ../Fltk/Opengl_Window.h \
+  ../Fltk/Colorbar_Window.h ../Common/GmshUI.h ../Fltk/Popup_Button.h \
   ../Fltk/SpherePosition_Widget.h ../Mesh/Field.h ../Geo/Geo.h \
   ../Common/GmshDefines.h ../Geo/gmshSurface.h ../Geo/Pair.h \
   ../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
diff --git a/Common/Message.cpp b/Common/Message.cpp
index f4862bda83700260f2b4a558c399c53309adaa27..29897ff1cc7c4da65fca40dcb557ffcd20fe1c3b 100644
--- a/Common/Message.cpp
+++ b/Common/Message.cpp
@@ -1,4 +1,4 @@
-// $Id: Message.cpp,v 1.4 2008-07-02 17:40:56 geuzaine Exp $
+// $Id: Message.cpp,v 1.5 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -23,6 +23,7 @@
 #include <stdio.h>
 #include <string.h>
 #include "Message.h"
+#include "Gmsh.h"
 #include "Options.h"
 #include "Context.h"
 #include "OS.h"
@@ -47,6 +48,7 @@ int Message::_progressMeterCurrent = 0;
 std::map<std::string, double> Message::_timers;
 int Message::_warningCount = 0;
 int Message::_errorCount = 0;
+GmshMessage *Message::_callback = 0;
 
 #if defined(HAVE_NO_VSNPRINTF)
 static int vsnprintf(char *str, size_t size, const char *fmt, va_list ap)
@@ -140,15 +142,17 @@ void Message::Fatal(const char *fmt, ...)
 {
   _errorCount++;
 
-  if(_verbosity < 1) return;
+  char str[1024];
   va_list args;
   va_start(args, fmt);
+  vsnprintf(str, sizeof(str), fmt, args);
+  va_end(args);
+
+  if(_callback) (*_callback)("Fatal", str);
 
 #if defined(HAVE_FLTK)
   if(WID){
     WID->check();
-    char str[1024];
-    vsnprintf(str, sizeof(str), fmt, args);
     std::string tmp = std::string("@C1@.") + "Fatal   : " + str;
     WID->add_message(tmp.c_str());
     WID->create_message_window();
@@ -159,16 +163,14 @@ void Message::Fatal(const char *fmt, ...)
 
   if(CTX.terminal){
     if(_commSize > 1)
-      fprintf(stderr, "Fatal   : [On processor %d] ", _commRank);
+      fprintf(stderr, "Fatal   : [On processor %d] %s\n", _commRank, str);
     else
-      fprintf(stderr, "Fatal   : ");
-    vfprintf(stderr, fmt, args);
-    fprintf(stderr, "\n");
+      fprintf(stderr, "Fatal   : %s\n", str);
     fflush(stderr);
   }
 
-  va_end(args);
-  Exit(1);
+  // only exit if a callback is not provided
+  if(!_callback) Exit(1);
 }
 
 void Message::Error(const char *fmt, ...)
@@ -176,15 +178,18 @@ void Message::Error(const char *fmt, ...)
   _errorCount++;
 
   if(_verbosity < 1) return;
+
+  char str[1024];
   va_list args;
   va_start(args, fmt);
+  vsnprintf(str, sizeof(str), fmt, args);
+  va_end(args);
+
+  if(_callback) (*_callback)("Error", str);
 
 #if defined(HAVE_FLTK)
   if(WID){
     WID->check();
-    char str[1024];
-    vsnprintf(str, sizeof(str), fmt, args);
-    va_end(args);
     std::string tmp = std::string("@C1@.") + "Error   : " + str;
     WID->add_message(tmp.c_str());
     WID->create_message_window();
@@ -193,15 +198,11 @@ void Message::Error(const char *fmt, ...)
 
   if(CTX.terminal){
     if(_commSize > 1) 
-      fprintf(stderr, "Error   : [On processor %d]", _commRank);
+      fprintf(stderr, "Error   : [On processor %d] %s\n", _commRank, str);
     else
-      fprintf(stderr, "Error   : ");
-    vfprintf(stderr, fmt, args);
-    fprintf(stderr, "\n");
+      fprintf(stderr, "Error   : %s\n", str);
     fflush(stderr);
   }
-
-  va_end(args);
 }
 
 void Message::Warning(const char *fmt, ...)
@@ -209,61 +210,62 @@ void Message::Warning(const char *fmt, ...)
   _warningCount++;
 
   if(_commRank || _verbosity < 2) return;
+
+  char str[1024];
   va_list args;
   va_start(args, fmt);
+  vsnprintf(str, sizeof(str), fmt, args);
+  va_end(args);
+
+  if(_callback) (*_callback)("Warning", str);
 
 #if defined(HAVE_FLTK)
   if(WID){
     WID->check();
-    char str[1024];
-    vsnprintf(str, sizeof(str), fmt, args);
     std::string tmp = std::string("@C1@.") + "Warning : " + str;
     WID->add_message(tmp.c_str());
   }
 #endif
 
   if(CTX.terminal){
-    fprintf(stderr, "Warning : ");
-    vfprintf(stderr, fmt, args);
-    fprintf(stderr, "\n");
+    fprintf(stderr, "Warning : %s\n", str);
     fflush(stderr);
   }
-
-  va_end(args);
 }
 
 void Message::Info(const char *fmt, ...)
 {
   if(_commRank || _verbosity < 3) return;
+
+  char str[1024];
   va_list args;
   va_start(args, fmt);
+  vsnprintf(str, sizeof(str), fmt, args);
+  va_end(args);
+
+  if(_callback) (*_callback)("Info", str);
 
 #if defined(HAVE_FLTK)
   if(WID){
     WID->check();
-    char str[1024];
-    vsnprintf(str, sizeof(str), fmt, args);
     std::string tmp = std::string("Info    : ") + str;
     WID->add_message(tmp.c_str());
   }
 #endif
 
   if(CTX.terminal){
-    fprintf(stdout, "Info    : ");
-    vfprintf(stdout, fmt, args);
-    fprintf(stdout, "\n");
+    fprintf(stdout, "Info    : %s\n", str);
     fflush(stdout);
   }
-
-  va_end(args);
 }
 
 void Message::Direct(const char *fmt, ...)
 {
   if(_commRank || _verbosity < 3) return;
+
+  char str[1024];
   va_list args;
   va_start(args, fmt);
-  char str[1024];
   vsnprintf(str, sizeof(str), fmt, args);
   va_end(args);
   Direct(3, str);
@@ -272,14 +274,18 @@ void Message::Direct(const char *fmt, ...)
 void Message::Direct(int level, const char *fmt, ...)
 {
   if(_commRank || _verbosity < level) return;
+
+  char str[1024];
   va_list args;
   va_start(args, fmt);
+  vsnprintf(str, sizeof(str), fmt, args);
+  va_end(args);
+
+  if(_callback) (*_callback)("Direct", str);
 
 #if defined(HAVE_FLTK)
   if(WID){
     WID->check();
-    char str[1024];
-    vsnprintf(str, sizeof(str), fmt, args);
     std::string tmp;
     if(level < 3)
       tmp = std::string("@C1@.") + str;
@@ -292,26 +298,27 @@ void Message::Direct(int level, const char *fmt, ...)
 #endif
   
   if(CTX.terminal){
-    vfprintf(stdout, fmt, args);
-    fprintf(stdout, "\n");
+    fprintf(stdout, "%s\n", str);
     fflush(stdout);
   }
-
-  va_end(args);
 }
 
 void Message::StatusBar(int num, bool log, const char *fmt, ...)
 {
   if(_commRank || _verbosity < 3) return;
   if(num < 1 || num > 3) return;
+
+  char str[1024];
   va_list args;
   va_start(args, fmt);
+  vsnprintf(str, sizeof(str), fmt, args);
+  va_end(args);
+
+  if(_callback && log) (*_callback)("Info", str);
 
 #if defined(HAVE_FLTK)
   if(WID){
     if(log) WID->check();
-    char str[1024];
-    vsnprintf(str, sizeof(str), fmt, args);
     WID->set_status(str, num - 1);
     if(log){
       std::string tmp = std::string("Info    : ") + str;
@@ -321,25 +328,25 @@ void Message::StatusBar(int num, bool log, const char *fmt, ...)
 #endif
 
   if(log && CTX.terminal){
-    fprintf(stdout, "Info    : ");
-    vfprintf(stdout, fmt, args);
-    fprintf(stdout, "\n");
+    fprintf(stdout, "Info    : %s\n", str);
     fflush(stdout);
   }
-
-  va_end(args);
 }
 
 void Message::Debug(const char *fmt, ...)
 {
   if(_verbosity < 99) return;
+
+  char str[1024];
   va_list args;
   va_start(args, fmt);
+  vsnprintf(str, sizeof(str), fmt, args);
+  va_end(args);
+
+  if(_callback) (*_callback)("Debug", str);
 
 #if defined(HAVE_FLTK)
   if(WID){
-    char str[1024];
-    vsnprintf(str, sizeof(str), fmt, args);
     std::string tmp = std::string("Debug   : ") + str;
     WID->add_message(tmp.c_str());
   }
@@ -347,15 +354,11 @@ void Message::Debug(const char *fmt, ...)
 
   if(CTX.terminal){
     if(_commSize > 1) 
-      fprintf(stdout, "Debug   : [On processor %d] ", _commRank);
+      fprintf(stdout, "Debug   : [On processor %d] %s\n", _commRank, str);
     else
-      fprintf(stdout, "Debug   : ");
-    vfprintf(stdout, fmt, args);
-    fprintf(stdout, "\n");
+      fprintf(stdout, "Debug   : %s\n", str);
     fflush(stdout);
   }
-
-  va_end(args);
 }
 
 void Message::ProgressMeter(int n, int N, const char *fmt, ...)
@@ -363,27 +366,29 @@ void Message::ProgressMeter(int n, int N, const char *fmt, ...)
   if(_commRank || _verbosity < 3) return;
 
   double percent = 100. * (double)n/(double)N;
+
   if(percent >= _progressMeterCurrent){
+    char str[1024], str2[1024];
     va_list args;
     va_start(args, fmt);
+    vsnprintf(str, sizeof(str), fmt, args);
+    va_end(args);
+
+    if(strlen(fmt)) strcat(str, " ");
+    sprintf(str2, "(%d %%)", _progressMeterCurrent);
+    strcat(str, str2);
+
 #if defined(HAVE_FLTK)
     if(WID){
-      char str[1024], str2[1024];
-      vsnprintf(str, sizeof(str), fmt, args);
-      if(strlen(fmt)) strcat(str, " ");
-      sprintf(str2, "(%d %%)", _progressMeterCurrent);
-      strcat(str, str2);
       WID->set_status(str, 1);
       WID->check();
     }
 #endif
     if(CTX.terminal){
-      vfprintf(stdout, fmt, args);
-      if(strlen(fmt)) fprintf(stdout, " ");
-      fprintf(stdout, "(%d %%)                     \r", _progressMeterCurrent);
+      fprintf(stdout, "%s                     \r", str);
       fflush(stdout);
     }
-    va_end(args);
+
     while(_progressMeterCurrent < percent)
       _progressMeterCurrent += _progressMeterStep;
   }
diff --git a/Common/Message.h b/Common/Message.h
index 7c519a4a53968f6b449866aac2a38a14d5cc27f6..95a8dcdb8d078835fb526428a5050218bbac0f38 100644
--- a/Common/Message.h
+++ b/Common/Message.h
@@ -25,6 +25,8 @@
 #include <string>
 #include <stdarg.h>
 
+class GmshMessage;
+
 // a class to manage messages
 class Message {
  private:
@@ -38,6 +40,8 @@ class Message {
   static std::map<std::string, double> _timers;
   // counters
   static int _warningCount, _errorCount;
+  // callback
+  static GmshMessage *_callback;
  public:
   Message() {}
   static void Init(int argc, char **argv);
@@ -46,6 +50,7 @@ class Message {
   static int GetCommSize(){ return _commSize; }
   static void SetCommRank(int val){ _commRank = val; }
   static void SetCommSize(int val){ _commSize = val; }
+  static void SetCallback(GmshMessage *callback){ _callback = callback; }
   static void Barrier();
   static int GetNumThreads();
   static int GetMaxThreads();
diff --git a/Common/SmoothData.cpp b/Common/SmoothData.cpp
index 0e575487be4f3be74e3bd19c8571b45196f9f4f6..eca379b201248d7f57dae62154fdf738d1b2a4bc 100644
--- a/Common/SmoothData.cpp
+++ b/Common/SmoothData.cpp
@@ -1,4 +1,4 @@
-// $Id: SmoothData.cpp,v 1.7 2008-06-05 11:52:49 samtech Exp $
+// $Id: SmoothData.cpp,v 1.8 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -67,9 +67,8 @@ void xyzv::update(int n, double *v)
     nbvals = n;
     nboccurences = 0;
   }
-  else if(nbvals != n) {
-    throw n;
-  }
+  else if(nbvals != n)
+    return; // error
   double x1 = (double)(nboccurences) / (double)(nboccurences + 1);
   double x2 = 1. / (double)(nboccurences + 1);
   for(int i = 0; i < nbvals; i++)
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 2600a17868991f96ea8f85f2ada3601245b08b01..f6e149753469c5e23a9a6a3693024737013b8b31 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdgeLoop.cpp,v 1.15 2008-06-05 11:52:49 samtech Exp $
+// $Id: GEdgeLoop.cpp,v 1.16 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -32,7 +32,7 @@
 void GEdgeSigned::print() const
 {
   Msg::Info("GEdgeSigned : Edge %d sign %d Ordered Vertices %d,%d",
-      ge->tag(), _sign, getBeginVertex()->tag(), getEndVertex()->tag());
+	    ge->tag(), _sign, getBeginVertex()->tag(), getEndVertex()->tag());
 }
 
 int countInList(std::list<GEdge*> &wire, GEdge *ge)
@@ -81,7 +81,7 @@ GEdgeSigned nextOne(GEdgeSigned *thisOne, std::list<GEdge*> &wire)
       GVertex *v2 = ge->getEndVertex();
       if(v1 == gv) return GEdgeSigned(1, ge);   
       if(v2 == gv) return GEdgeSigned(-1, ge);   
-      throw;
+      Msg::Error("Something wrong in edge loop");
     }
     ++it;
   }
@@ -97,7 +97,7 @@ GEdgeSigned nextOne(GEdgeSigned *thisOne, std::list<GEdge*> &wire)
       GVertex *v2 = ge->getEndVertex();
       if(v1 == gv) return GEdgeSigned(1, ge);   
       if(v2 == gv) return GEdgeSigned(-1, ge);
-      throw;
+      Msg::Error("Something wrong in edge loop");
     }   
     ++it;
   }
@@ -124,16 +124,15 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
 
   GEdgeSigned *prevOne = 0;
 
-  //  Msg::Info("Building a wire");
   GEdgeSigned ges(0,0);
   while(wire.size()){
     ges = nextOne(prevOne, wire);
     if(ges.getSign() == 0){ // oops
-      Msg::Warning("Something wrong in edge loop?");
+      Msg::Error("Something wrong in edge loop");
       break;
     }
     prevOne = &ges;
-    //    ges.print();
+    // ges.print();
     loop.push_back(ges);
   }
 }
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index f557f89ac968c66f9aa84bccf5bec49361c36b1e..d03bfec6e7e59763532461d9f21a13a8457af25c 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -161,22 +161,22 @@ class GEntity {
   void deleteVertexArrays();
 
   // Spatial dimension of the entity
-  virtual int dim() const { throw; }
+  virtual int dim() const { return -1; }
 
   // Regions that bound this entity or that this entity bounds.
-  virtual std::list<GRegion*> regions() const { throw; }
+  virtual std::list<GRegion*> regions() const { return std::list<GRegion*>(); }
 
   // Faces that bound this entity or that this entity bounds.
-  virtual std::list<GFace*> faces() const { throw; }
+  virtual std::list<GFace*> faces() const { return std::list<GFace*>(); }
 
   // Edges that bound this entity or that this entity bounds.
-  virtual std::list<GEdge*> edges() const { throw; }
+  virtual std::list<GEdge*> edges() const { return std::list<GEdge*>(); }
 
   // Vertices that bound this entity.
-  virtual std::list<GVertex*> vertices() const { throw; }
+  virtual std::list<GVertex*> vertices() const { return std::list<GVertex*>(); }
 
   /// Underlying geometric representation of this entity.
-  virtual GeomType geomType() const { throw; }
+  virtual GeomType geomType() const { return Unknown; }
 
   // True if parametric space is continuous in the "dim" direction.
   virtual bool continuous(int dim) const { return true; }
@@ -188,19 +188,19 @@ class GEntity {
   virtual bool degenerate(int dim) const { return false; }
 
   // Parametric bounds of the entity in the "i" direction.
-  virtual Range<double> parBounds(int i) const { throw; }
+  virtual Range<double> parBounds(int i) const { return Range<double>(0., 0.); }
 
   // Modeler tolerance for the entity.
   virtual double tolerance() const { return 1.e-14; }
 
   // True if the entity contains the given point to within tolerance.
-  virtual int containsPoint(const SPoint3 &pt) const { throw; }
+  virtual bool containsPoint(const SPoint3 &pt) const { return false; }
 
   // Get the native type of the particular representation
   virtual ModelType getNativeType() const { return UnknownModel; }
 
   // Get the native pointer of the particular representation
-  virtual void * getNativePtr() const { throw; }
+  virtual void * getNativePtr() const { return 0; }
 
   // The model owning this entity.
   GModel *model() const { return _model; }
@@ -210,7 +210,7 @@ class GEntity {
   void setTag(int tag) { _tag = tag; }
 
   // The bounding box
-  virtual SBoundingBox3d bounds() const { throw; }
+  virtual SBoundingBox3d bounds() const { return SBoundingBox3d(); }
 
   // Get the visibility flag
   virtual char getVisibility();
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 30241dc5c7c44377892598652c893b3f003275b1..c136917ffc97defe66790f5c040161460a6c0f73 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.92 2008-07-02 17:40:56 geuzaine Exp $
+// $Id: GModel.cpp,v 1.93 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -208,7 +208,10 @@ void GModel::snapVertices()
       else if ((*it)->getEndVertex() == *vit){
         t = parb.high();
       }
-      else throw;
+      else{
+	Msg::Error("Weird vertex: impossible to snap");
+	break;
+      }
       GPoint gp = (*it)->point(t);
       double d = sqrt((gp.x() - (*vit)->x()) * (gp.x() - (*vit)->x()) +
                       (gp.y() - (*vit)->y()) * (gp.y() - (*vit)->y()) +
@@ -216,7 +219,7 @@ void GModel::snapVertices()
       if (d > tol){
         (*vit)->setPosition(gp);
         Msg::Warning("Geom Vertex %d Corrupted (%12.5E)... Snap performed",
-            (*vit)->tag(), d);
+		     (*vit)->tag(), d);
       }
     }
     vit++;
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index bc24ddd68067a4a5ae665a7b68d863a14115edea..53cd8a6562d01f01e63b71203e2cb6317fcf091b 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_OCC.cpp,v 1.39 2008-07-02 17:40:56 geuzaine Exp $
+// $Id: GModelIO_OCC.cpp,v 1.40 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -526,18 +526,17 @@ void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool, const BooleanOperato
             if (aValueC.ShapeType() != TopAbs_SOLID) isOnlySolids = false;
           }
           if (isOnlySolids)
-            throw;
-            //      theNewShape = GEOMImpl_GlueDriver::GlueFaces(C, Precision::Confusion());
-          else
-            theNewShape = C;
+	    Msg::Error("Face gluing not implemented");
+	  //  theNewShape = GEOMImpl_GlueDriver::GlueFaces(C, Precision::Confusion());
+          //else
+	      theNewShape = C;
         }       
       }
       break;
     case OCC_Internals::Cut :
-      {
-      }
     default :
-      throw;
+      Msg::Error("Requested boolean operation not implemented");
+      break;
     }
   }
 }
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 8c663ea2812720b4f6cf0b0ab2259d6d69e9046c..823ae6e93dfbec18fdff16bbf1445abb4317fbd0 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.39 2008-06-27 17:34:19 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.40 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -101,14 +101,12 @@ static Vertex InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
 
 // interpolation in the parametric space !
 SPoint2 InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
-                              int derivee, double t1, double t2, gmshSurface *s)
+			       double t1, double t2, gmshSurface *s)
 {
   Vertex V;
   int i, j;
   double T[4];
 
-  if(derivee) throw;
-
   T[3] = 1.;
   T[2] = t;
   T[1] = t * t;
@@ -148,7 +146,7 @@ static Vertex InterpolateUBS(Curve *Curve, double u, int derivee)
   }
 
   if(Curve->geometry){
-    SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, 0, t1, t2, Curve->geometry);
+    SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, t1, t2, Curve->geometry);
     SPoint3 pt = Curve->geometry->point(pp);
     Vertex V;
     V.Pos.X = pt.x();
@@ -366,7 +364,7 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
       List_Read(c->Control_Points, i + 2, &v[3]);
     }
     if(c->geometry){
-      SPoint2 pp = InterpolateCubicSpline(v, t, c->mat, 0, t1, t2,c->geometry);
+      SPoint2 pp = InterpolateCubicSpline(v, t, c->mat, t1, t2,c->geometry);
       SPoint3 pt = c->geometry->point(pp);
       V.Pos.X = pt.x();
       V.Pos.Y = pt.y();
diff --git a/Geo/GeoInterpolation.h b/Geo/GeoInterpolation.h
index aa00d52797f7ccb3b1fa2d5fbec24fa68dae9cc3..4978e69804d9c8f18165b86c657f2ba2984c2679 100644
--- a/Geo/GeoInterpolation.h
+++ b/Geo/GeoInterpolation.h
@@ -26,5 +26,5 @@ bool iSRuledSurfaceASphere(Surface *s, SPoint3 &center, double &radius);
 Vertex InterpolateCurve(Curve *Curve, double u, int derivee);
 Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v);
 SPoint2 InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
-                               int derivee, double t1, double t2, gmshSurface *s);
+                               double t1, double t2, gmshSurface *s);
 #endif
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 92c620ddb3ff9fe3e1a2c748d3cf1b2787e00a16..4e85bfe43f4d8bd5b5c0b8326b31686940dd7b46 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.76 2008-06-27 08:33:36 koen Exp $
+// $Id: MElement.cpp,v 1.77 2008-07-03 17:06:01 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -105,6 +105,11 @@ double MElement::rhoShapeMeasure()
     return 0.;
 }
 
+void MElement::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+{
+  Msg::Error("No integration points defined for this type of element");
+}
+
 double MTriangle::gammaShapeMeasure()
 {
 #if defined(HAVE_GMSH_EMBEDDED)
@@ -647,8 +652,8 @@ void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double ww, dou
     case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
     case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu, vv, ww, grads); break;
     case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu, vv, ww, grads); break;
-    case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv,ww, grads); break;
-    default: throw;
+    case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv, ww, grads); break;
+    default: Msg::Error("Order %d triangle jac not implemented", ord); break;
     }
   }
   else{
@@ -658,7 +663,7 @@ void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double ww, dou
     case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu, vv, ww,grads); break;
     case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu, vv, ww,grads); break;
     case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu, vv, ww,grads); break;
-    default: throw;
+    default: Msg::Error("Order %d triangle jac not implemented", ord); break;
     }
   }
   j[0][0] = 0 ; for(int i = 0; i < 3; i++) j[0][0] += grads [i][0] * _v[i]->x();
@@ -692,7 +697,7 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, double ww, SPo
     case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, ww,sf); break;
     case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, ww,sf); break;
     case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, ww,sf); break;
-    default: throw;
+    default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
     }
   }
   else{
@@ -702,7 +707,7 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, double ww, SPo
     case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, ww,sf); break;
     case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, ww,sf); break;
     case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, ww,sf); break;
-    default: throw;
+    default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
     }
   }
   
@@ -728,7 +733,7 @@ void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, double ww,S
   case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww,sf); break;
   case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww,sf); break;
   case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww,sf); break;
-  default : throw;
+  default: Msg::Error("Order %d tetrahedron pnt not implemented", ord); break;
   }
     
   double x = 0 ; for(int i = 0; i < 4; i++) x += sf[i] * _v[i]->x();
@@ -745,6 +750,7 @@ void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, double ww,S
   p = SPoint3(x,y,z);
 #endif
 }
+
 void MTetrahedron::pnt(int ord, std::vector<MVertex *> & vs, double uu, double vv, double ww,SPoint3 &p)
 {
 #if !defined(HAVE_GMSH_EMBEDDED)
@@ -755,7 +761,7 @@ void MTetrahedron::pnt(int ord, std::vector<MVertex *> & vs, double uu, double v
   case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww,sf); break;
   case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww,sf); break;
   case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww,sf); break;
-  default : throw;
+  default: Msg::Error("Order %d tetrahedron pnt not implemented", ord); break;
   }
     
   double x = 0 ; for(int i = 0; i < 4; i++) x += sf[i] * _v[i]->x();
@@ -791,7 +797,7 @@ void MTetrahedron::jac(int ord, MVertex *vs[], double uu, double vv, double ww,
   case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, ww, grads); break;
   case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, ww, grads); break;
   case 5: gmshFunctionSpaces::find(MSH_TET_56).df(uu, vv, ww, grads); break;
-  default: throw;
+  default: Msg::Error("Order %d tetrahedron jac not implemented", ord); break;
   }
  
   j[0][0] = 0 ; for(int i = 0; i < 4; i++) j[0][0] += grads [i][0] * _v[i]->x();
@@ -838,7 +844,7 @@ void MTetrahedron::jac(int ord, std::vector<MVertex *>& vs, double uu, double vv
   case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, ww, grads); break;
   case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, ww, grads); break;
   case 5: gmshFunctionSpaces::find(MSH_TET_56).df(uu, vv, ww, grads); break;
-  default: throw;
+  default: Msg::Error("Order %d tetrahedron jac not implemented", ord); break;
   }
  
   j[0][0] = 0 ; for(int i = 0; i < 4; i++) j[0][0] += grads [i][0] * _v[i]->x();
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 65702bca0d2e920adf66e4b60e78d13fe11e8bb5..46bde832b2498c2ade1996e841b58d19eeafd385 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -189,7 +189,7 @@ class MElement
                        int stride=3);
   double interpolateDiv(double val[], double u, double v, double w, int stride=3);
   // integration routine 
-  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const { throw; }
+  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
   
   // IO routines
   virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, 
@@ -260,12 +260,9 @@ class MLine : public MElement {
     _getEdgeRep(_v[0], _v[1], x, y, z, n);
   }
   virtual int getNumFaces(){ return 0; }
-  virtual MFace getFace(int num){ throw; }
+  virtual MFace getFace(int num){ return MFace(); }
   virtual int getNumFacesRep(){ return 0; }
-  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
-  { 
-    throw; 
-  }
+  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){}
   virtual int getTypeForMSH(){ return MSH_LIN_2; }
   virtual int getTypeForUNV(){ return 21; } // linear beam
   virtual int getTypeForVTK(){ return 3; }
@@ -297,7 +294,7 @@ class MLine : public MElement {
       return false;
     return true;
   }
-  virtual void getIntegrationPoints ( int pOrder , int *npts, IntPt **pts) const;
+  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
 };
 
 class MLine3 : public MLine {
@@ -407,7 +404,7 @@ class MTriangle : public MElement {
     if(_v[0] != v1 && _v[0] != v2) return _v[0];
     if(_v[1] != v1 && _v[1] != v2) return _v[1];
     if(_v[2] != v1 && _v[2] != v2) return _v[2];
-    throw;
+    return 0;
   }
   virtual int getNumEdges(){ return 3; }
   virtual MEdge getEdge(int num)
@@ -478,7 +475,7 @@ class MTriangle : public MElement {
       return false; 
     return true;
   }
-  virtual void getIntegrationPoints ( int pOrder , int *npts, IntPt **pts) const;
+  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
 };
 
 class MTriangle6 : public MTriangle {
@@ -591,7 +588,7 @@ class MTriangleN : public MTriangle {
     if(_order == 4 && _vs.size() == 12) return 3;
     if(_order == 5 && _vs.size() == 12) return 0;
     if(_order == 5 && _vs.size() == 18) return 6;
-    throw;
+    return 0;
   }
   virtual int getNumEdgeVertices(){ return _order - 1; }
   virtual int getNumEdgesRep();
@@ -712,7 +709,7 @@ class MQuadrangle : public MElement {
       return false;
     return true;
   }
-  virtual void getIntegrationPoints ( int pOrder , int *npts, IntPt **pts) const;
+  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
 };
 
 class MQuadrangle8 : public MQuadrangle {
@@ -957,7 +954,7 @@ class MTetrahedron : public MElement {
       return false;
     return true;
   }
-  virtual void getIntegrationPoints ( int pOrder , int *npts, IntPt **pts) const;
+  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
 };
 
 class MTetrahedron10 : public MTetrahedron {
@@ -1075,8 +1072,7 @@ class MTetrahedronN : public MTetrahedron {
     case MSH_TET_20 : return 1;
     case MSH_TET_35 : return 3;
     case MSH_TET_56 : return 6;
-    default : 
-      throw;
+    default : return 0;
     }    
   }
   virtual int getNumEdgeVertices(){ return _order - 1; }
@@ -1086,7 +1082,7 @@ class MTetrahedronN : public MTetrahedron {
     if(_order == 3 && _vs.size() + 4 == 20) return MSH_TET_20; 
     if(_order == 4 && _vs.size() + 4 == 35) return MSH_TET_35; 
     if(_order == 5 && _vs.size() + 4 == 56) return MSH_TET_56; 
-    throw;
+    return 0;
   }
   virtual void revert() 
   {
@@ -1251,7 +1247,7 @@ class MHexahedron : public MElement {
       return false;
     return true;
   }
-  virtual void getIntegrationPoints ( int pOrder , int *npts, IntPt **pts) const;
+  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
 };
 
 class MHexahedron20 : public MHexahedron {
diff --git a/Geo/MFace.cpp b/Geo/MFace.cpp
index 795c38608f3e350fec2ee5d8a8ebb2e8ccfb64e8..be8a2efca7ebc978a36bde2ff10da6642c97cfe9 100644
--- a/Geo/MFace.cpp
+++ b/Geo/MFace.cpp
@@ -28,6 +28,14 @@
 
 extern Context_T CTX;
 
+MFace::MFace() 
+{ 
+  for(int i = 0; i < 4; i++){
+    _v[i] = 0; 
+    _si[i] = 0;
+  }
+}
+
 MFace::MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3) 
 {
   if(CTX.mesh.reverse_all_normals){
diff --git a/Geo/MFace.h b/Geo/MFace.h
index 2a0146ca35794982f59e0e7e511b14ab823c195a..e86ce29b4538b5b88292d221ebf4d1b77d11ea61 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -32,6 +32,7 @@ class MFace {
   char _si[4]; // sorted indices
 
  public:
+  MFace();
   MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3=0);
   inline int getNumVertices() const { return _v[3] ? 4 : 3; }
   inline MVertex *getVertex(const int i) const { return _v[i]; }
@@ -100,7 +101,6 @@ class MFace {
         p[2] += v->z() * ff[i] * .25;
       } 
     }
-    else throw;
     return p;
   }
 };
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 3375b0627b5e73ffdea8dfbd5dc5adb1eb10bea5..46f66e662b9b82d959258cf313f34ca0558de8dd 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -141,13 +141,11 @@ class MEdgeVertex : public MVertex{
   virtual ~MEdgeVertex(){}
   virtual bool getParameter(int i, double &par) const 
   { 
-    if(i) throw; 
     par = _u; 
     return true; 
   }
   virtual bool setParameter(int i, double par)
   { 
-    if(i) throw; 
     _u = par; 
     return true; 
   }
diff --git a/Geo/Makefile b/Geo/Makefile
index ddee1cbd9a2b9dac7c5d3e4c86db1871e28e45b5..bca6d43259ed27153bb2f72c05acfd6fd0d0d19a 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.215 2008-07-01 15:11:38 geuzaine Exp $
+# $Id: Makefile,v 1.216 2008-07-03 17:06:02 geuzaine Exp $
 #
 # Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 #
@@ -161,14 +161,14 @@ discreteEdge.o: discreteEdge.cpp discreteEdge.h GModel.h GVertex.h \
   ../Common/GmshDefines.h gmshSurface.h ../Numeric/Numeric.h \
   ../Numeric/NumericEmbedded.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
-  ExtrudeParams.h ../Common/SmoothData.h
+  ExtrudeParams.h ../Common/SmoothData.h ../Common/Message.h
 discreteFace.o: discreteFace.cpp discreteFace.h GModel.h GVertex.h \
   GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h ../Numeric/Numeric.h \
   ../Numeric/NumericEmbedded.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
-  ExtrudeParams.h ../Common/SmoothData.h
+  ExtrudeParams.h ../Common/SmoothData.h ../Common/Message.h
 discreteRegion.o: discreteRegion.cpp discreteRegion.h GModel.h GVertex.h \
   GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h Geo.h \
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 30016b76f1a37779036b7382a261a4c1adb8f30c..798dd975236ee54844ff489641625002b9e43a4d 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.38 2008-05-04 08:31:13 geuzaine Exp $
+// $Id: OCCEdge.cpp,v 1.39 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -77,7 +77,7 @@ SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar, int dir) const
   
   if(c2d.IsNull()){
     Msg::Fatal("Reparam on face failed: curve %d is not on surface %d",
-        tag(), face->tag());
+	       tag(), face->tag());
   }
 
   double u, v;
@@ -93,11 +93,11 @@ SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar, int dir) const
   if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-4 * CTX.lc){
     // return reparamOnFace(face, epar,-1);      
     Msg::Warning("Reparam on face partially failed for curve %d surface %d at point %g",
-        tag(), face->tag(), epar);
+		 tag(), face->tag(), epar);
     Msg::Warning("On the face %d local (%g %g) global (%g %g %g)",
-        face->tag(), u, v, p2.x(), p2.y(), p2.z());
+		 face->tag(), u, v, p2.x(), p2.y(), p2.z());
     Msg::Warning("On the edge %d local (%g) global (%g %g %g)",
-        tag(), epar, p1.x(), p1.y(), p1.z());
+		 tag(), epar, p1.x(), p1.y(), p1.z());
     // GPoint ppp = face->closestPoint(SPoint3(p1.x(), p1.y(), p1.z()));
     // return SPoint2(ppp.u(), ppp.v());
   }
@@ -145,7 +145,8 @@ SVector3 OCCEdge::firstDer(double par) const
 
 double OCCEdge::parFromPoint(const SPoint3 &pt) const
 {
-  throw;
+  Msg::Error("parFromPoint not implemented for OCCEdge");
+  return 0.;
 }
 
 GEntity::GeomType OCCEdge::geomType() const
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index e95367c9dc294ccc4a634d2932549ea9b507cec1..bab19e83c123c4cf71620e10462cfa3ad4d47235 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.42 2008-06-10 08:37:34 remacle Exp $
+// $Id: OCCFace.cpp,v 1.43 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -46,28 +46,29 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
   TopExp_Explorer exp0, exp01, exp1, exp2, exp3;
   for(exp2.Init(s, TopAbs_WIRE); exp2.More(); exp2.Next()){
     TopoDS_Shape wire = exp2.Current();
-    Msg::Debug("OCC Face %d - New Wire",num);
+    Msg::Debug("OCC Face %d - New Wire", num);
     std::list<GEdge*> l_wire;
-    std::list<int> l_oris;
     for(exp3.Init(wire, TopAbs_EDGE); exp3.More(); exp3.Next()){          
       TopoDS_Edge edge = TopoDS::Edge(exp3.Current());
       int index = emap.FindIndex(edge);
       GEdge *e = m->getEdgeByTag(index);
-      if(!e) throw;
-      l_wire.push_back(e);
-      l_oris.push_back(edge.Orientation());
-      Msg::Debug("Edge %d ori %d", e->tag(), edge.Orientation());
-      e->addFace(this);
-      if(!e->is3D()){
-        OCCEdge *occe = (OCCEdge*)e;
-        occe->setTrimmed(this);
+      if(e){
+	l_wire.push_back(e);
+	Msg::Debug("Edge %d ori %d", e->tag(), edge.Orientation());
+	e->addFace(this);
+	if(!e->is3D()){
+	  OCCEdge *occe = (OCCEdge*)e;
+	  occe->setTrimmed(this);
+	}
+      }
+      else{
+	Msg::Error("Unknown edge %d in face %d", index, num);
       }
     }      
     
-    //    GEdgeLoop el(l_wire,l_oris);
-    GEdgeLoop el(l_wire);        
+    GEdgeLoop el(l_wire);
 
-    for(GEdgeLoop::citer it = el.begin() ; it != el.end() ; ++it){
+    for(GEdgeLoop::citer it = el.begin(); it != el.end(); ++it){
       l_edges.push_back(it->ge);
       l_dirs.push_back(it->_sign);
       if (el.count() == 2){
@@ -214,7 +215,7 @@ double OCCFace::curvature (const SPoint2 &param) const
   return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
 }
 
-int OCCFace::containsPoint(const SPoint3 &pt) const
+bool OCCFace::containsPoint(const SPoint3 &pt) const
 { 
   if(geomType() == Plane){
     gp_Pln pl = Handle(Geom_Plane)::DownCast(occface)->Pln();
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index b3d2f3f5a597b01d8e3f900718d23088363c22a4..25662ebc096a4d192210722145aede3758fba4e0 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -41,7 +41,7 @@ class OCCFace : public GFace {
   Range<double> parBounds(int i) const; 
   virtual GPoint point(double par1, double par2) const; 
   virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const; 
-  virtual int containsPoint(const SPoint3 &pt) const;  
+  virtual bool containsPoint(const SPoint3 &pt) const;  
   virtual SVector3 normal(const SPoint2 &param) const; 
   virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const; 
   virtual GEntity::GeomType geomType() const; 
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index 35db2c91b01a8a73e444d841cfc126d3f631e79b..722cb1564804e3a590f81c03fc5d1e6b2925b3f4 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCRegion.cpp,v 1.12 2008-05-04 08:31:13 geuzaine Exp $
+// $Id: OCCRegion.cpp,v 1.13 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -39,9 +39,12 @@ OCCRegion::OCCRegion(GModel *m, TopoDS_Solid _s, int num, TopTools_IndexedMapOfS
       TopoDS_Face face = TopoDS::Face(exp3.Current());
       int index = fmap.FindIndex(face);
       GFace *f = m->getFaceByTag(index);
-      if(!f) throw;
-      l_faces.push_back(f);
-      f->addRegion(this);
+      if(f){
+	l_faces.push_back(f);
+	f->addRegion(this);
+      }
+      else
+	Msg::Error("Unknown face %d in region %d", index, num);
     }      
   }
   Msg::Info("OCC Region %d with %d edges", num, l_faces.size());
diff --git a/Geo/Pair.h b/Geo/Pair.h
index f31acea76a2ea1729248e0dccf618a6d587b4800..d98bdd204cd07c18675c1dd1e3e312bd1fd97c0e 100644
--- a/Geo/Pair.h
+++ b/Geo/Pair.h
@@ -28,7 +28,7 @@ class Pair{
   R Right;
  public:
   Pair() {}
-  Pair( const L& left, const R& right) : Left(left), Right(right) {}
+  Pair(const L& left, const R& right) : Left(left), Right(right) {}
   L left() const { return Left; }
   void left(const L& left) { Left = left; }
   R right() const { return Right; }
diff --git a/Geo/Range.h b/Geo/Range.h
index 3abbb0a09f8fd2afead7acb76f1911d2ffefc242..4bfe6ad473d86865ec67a6d5a1984826a9fbd2c8 100644
--- a/Geo/Range.h
+++ b/Geo/Range.h
@@ -28,7 +28,7 @@ class Range {
   T High;
  public:
   Range() {}
-  Range( const T& low, const T& high) : Low(low), High(high) {}
+  Range(const T& low, const T& high) : Low(low), High(high) {}
   T low() const { return Low; }
   void low(const T& low) { Low = low; }
   T high() const { return High; }
diff --git a/Geo/SPoint2.h b/Geo/SPoint2.h
index e89fe038e5c978f19ec8aac98d56bb16c89feec7..0e2b352f3382e53ba49518e81e76db4715aeaa60 100644
--- a/Geo/SPoint2.h
+++ b/Geo/SPoint2.h
@@ -27,9 +27,9 @@ class SPoint2 {
  protected:
   double P[2];
  public:
-  SPoint2(double x=0.0, double y=0.0) {P[0] = x; P[1] = y;}
-  SPoint2(double *p) {P[0] = p[0]; P[1] = p[1];}
-  SPoint2(const SPoint2 &pt) {P[0] = pt.P[0]; P[1] = pt.P[1]; }
+  SPoint2(double x=0.0, double y=0.0) { P[0] = x; P[1] = y; }
+  SPoint2(double *p) { P[0] = p[0]; P[1] = p[1]; }
+  SPoint2(const SPoint2 &pt) { P[0] = pt.P[0]; P[1] = pt.P[1]; }
   virtual ~SPoint2() {}
   void setPosition(double xx, double yy);
   void getPosition(double *xx, double *yy) const;
diff --git a/Geo/SPoint3.h b/Geo/SPoint3.h
index 2d05ae76d0bd2253e797baff8eb1fab24904e1d8..e80ad2ecadd503912f3791b69e9889dce8953bbc 100644
--- a/Geo/SPoint3.h
+++ b/Geo/SPoint3.h
@@ -27,7 +27,7 @@ class SPoint3 {
  protected:
   double P[3];
  public:
-  SPoint3() {}
+  SPoint3() { P[0] = P[1] = P[2] = 0.; }
   SPoint3(double x, double y, double z) { P[0] = x; P[1] = y; P[2] = z; }
   SPoint3(const double *p) { P[0] = p[0]; P[1] = p[1]; P[2] = p[2]; }
   SPoint3(const SPoint3 &pt) { P[0] = pt.P[0]; P[1] = pt.P[1]; P[2] = pt.P[2]; }
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index 5e473872777d9d29cee71347ea53112d7795b800..48ba5255e6054b1fc70510179f4c8ec6a89ce2db 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -24,6 +24,8 @@
 
 // concrete class for vector of size 3
 class SVector3 {
+ protected:
+  SPoint3 P;
  public:
   SVector3() {}
   // Construct from 2 SPoints, vector from p1 to p2
@@ -71,8 +73,6 @@ class SVector3 {
     return *this;
   }
   operator double *() { return P; }
- protected:
-  SPoint3 P;
 };
 
 inline double dot(const SVector3 &a, const SVector3 &b)
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 369a03a3a4cf72d1218031642d8b1f0613f2005d..7e4e64bf45124d51f488f436f29be5adc3defbd0 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: discreteEdge.cpp,v 1.1 2008-05-21 10:50:37 geuzaine Exp $
+// $Id: discreteEdge.cpp,v 1.2 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -21,8 +21,11 @@
 
 #include "discreteEdge.h"
 
-#if !defined(HAVE_GMSH_EMBEDDED)
-#include "Geo.h"
+#if defined(HAVE_GMSH_EMBEDDED)
+#  include "GmshEmbedded.h"
+#else
+#  include "Geo.h"
+#  include "Message.h"
 #endif
 
 discreteEdge::discreteEdge(GModel *model, int num) : GEdge(model, num, 0, 0) 
@@ -33,3 +36,21 @@ discreteEdge::discreteEdge(GModel *model, int num) : GEdge(model, num, 0, 0)
   CreateReversedCurve(c);
 #endif
 }
+
+GPoint discreteEdge::point(double p) const 
+{
+  Msg::Error("Cannot evaluate point on discrete edge");
+  return GPoint();
+}
+
+SVector3 discreteEdge::firstDer(double par) const 
+{
+  Msg::Error("Cannot evaluate derivative on discrete edge");
+  return SVector3();
+}
+
+double discreteEdge::parFromPoint(const SPoint3 &pt) const 
+{
+  Msg::Error("Cannot compute parametric coordinate on discrete edge");
+  return 0.;
+}
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index c8d86bbaeca22066386f1542255e828cf99e3b56..9831acc497b0ed5a449a66bf316b287b64f6b333 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -28,9 +28,9 @@ class discreteEdge : public GEdge {
   discreteEdge(GModel *model, int num);
   virtual ~discreteEdge() {}
   virtual GeomType geomType() const { return DiscreteCurve; }
-  virtual GPoint point(double p) const { throw; }
-  virtual SVector3 firstDer(double par) const { throw; }
-  virtual double parFromPoint(const SPoint3 &pt) const { throw; }
+  virtual GPoint point(double p) const;
+  virtual SVector3 firstDer(double par) const;
+  virtual double parFromPoint(const SPoint3 &pt) const;
 };
 
 #endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index b0f9e3a5abd1c279da8b0b275d5f5601b74503e7..e9e925fb68e0543978a7e99b457a3ef4ba286c3d 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -1,4 +1,4 @@
-// $Id: discreteFace.cpp,v 1.2 2008-05-20 19:25:33 geuzaine Exp $
+// $Id: discreteFace.cpp,v 1.3 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -21,8 +21,11 @@
 
 #include "discreteFace.h"
 
-#if !defined(HAVE_GMSH_EMBEDDED)
-#include "Geo.h"
+#if defined(HAVE_GMSH_EMBEDDED)
+#  include "GmshEmbedded.h"
+#else
+#  include "Geo.h"
+#  include "Message.h"
 #endif
 
 discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
@@ -33,3 +36,27 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
 #endif
   meshStatistics.status = GFace::DONE;    
 }
+
+GPoint discreteFace::point(double par1, double par2) const 
+{
+  Msg::Error("Cannot evaluate point on discrete face");
+  return GPoint();
+}
+
+SPoint2 discreteFace::parFromPoint(const SPoint3 &p) const
+{
+  Msg::Error("Cannot compute parametric coordinates on discrete face");
+  return SPoint2();
+}
+
+SVector3 discreteFace::normal(const SPoint2 &param) const
+{
+  Msg::Error("Cannot evaluate normal on discrete face");
+  return SVector3();
+}
+
+Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 &param) const
+{
+  Msg::Error("Cannot evaluate derivative on discrete face");
+  return Pair<SVector3, SVector3>();
+}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 47a7de6d052ed353d21ac70a5f40faa499af6f27..3de48fdc29be8a12bb92ddc694491b80e7e10db7 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -27,12 +27,11 @@ class discreteFace : public GFace {
  public:
   discreteFace(GModel *model, int num);
   virtual ~discreteFace() {}
-  virtual GPoint point(double par1, double par2) const { throw; }
-  virtual SPoint2 parFromPoint(const SPoint3 &p) const { throw; }
-  virtual SVector3 normal(const SPoint2 &param) const { throw; }
+  virtual GPoint point(double par1, double par2) const;
+  virtual SPoint2 parFromPoint(const SPoint3 &p) const;
+  virtual SVector3 normal(const SPoint2 &param) const;
   virtual GEntity::GeomType geomType() const { return DiscreteSurface; }
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const { throw; }
+  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const;
 };
 
 #endif
-
diff --git a/Geo/fourierFace.cpp b/Geo/fourierFace.cpp
index 9a4081f696202b5e7d6be3d079840c807f41184f..61da737eb62c7ce89064520a70dea552415300a5 100644
--- a/Geo/fourierFace.cpp
+++ b/Geo/fourierFace.cpp
@@ -59,4 +59,10 @@ GEntity::GeomType fourierFace::geomType() const
   return  GEntity::ParametricSurface;
 }
 
+Pair<SVector3, SVector3> fourierFace::firstDer(const SPoint2 &param) const
+{
+  Msg::Error("First derivative not implemented for fourier face");
+  return Pair<SVector3, SVector3>();
+}
+
 #endif
diff --git a/Geo/fourierFace.h b/Geo/fourierFace.h
index 538f7e55205c6196b0e96688620012b90368b200..8580a70c8530e7548e05ab8e74842ab33709f142 100644
--- a/Geo/fourierFace.h
+++ b/Geo/fourierFace.h
@@ -22,7 +22,7 @@ class fourierFace : public GFace {
   virtual int containsParam(const SPoint2 &pt) const; 
   virtual SVector3 normal(const SPoint2 &param) const; 
   virtual GEntity::GeomType geomType() const;
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const { throw; }
+  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const;
   ModelType getNativeType() const { return FourierModel; }
   void * getNativePtr() const { return face; } 
 };
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index a91b9605fece7b6333da64e4234ab5dd6bda8110..15b4142c99498491b1027695095536012d4680a0 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshEdge.cpp,v 1.49 2008-05-04 08:31:14 geuzaine Exp $
+// $Id: gmshEdge.cpp,v 1.50 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -168,7 +168,7 @@ SPoint2 gmshEdge::reparamOnFace(GFace *face, double epar,int dir) const
             k = periodic ? k - NbControlPoints + 1: NbControlPoints - 1;
           List_Read(c->Control_Points, k, &v[j]);
         }
-        return InterpolateCubicSpline(v, t, c->mat, 0, t1, t2, c->geometry);
+        return InterpolateCubicSpline(v, t, c->mat, t1, t2, c->geometry);
       }
     case MSH_SEGM_SPLN :
       {
@@ -209,10 +209,11 @@ SPoint2 gmshEdge::reparamOnFace(GFace *face, double epar,int dir) const
         else{
           List_Read(c->Control_Points, i + 2, &v[3]);
         }
-        return InterpolateCubicSpline(v, t, c->mat, 0, t1, t2, c->geometry);
+        return InterpolateCubicSpline(v, t, c->mat, t1, t2, c->geometry);
       }
     default:
-      throw;
+      Msg::Error("Unknown edge type in reparamOnFace");
+      return SPoint2();
     }
   }
   
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 6e22aa878528b7cfcc495235fbb652df82cd1e15..7ec6a112a6c593fbd2754f00f24e66dfcbee71f9 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.64 2008-07-01 14:24:07 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.65 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -268,7 +268,7 @@ SPoint2 gmshFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
   }
 }
 
-int gmshFace::containsPoint(const SPoint3 &pt) const
+bool gmshFace::containsPoint(const SPoint3 &pt) const
 { 
   if(geomType() == Plane){
     // OK to use the normal from the mean plane here: we compensate
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index 073caa7cbacb53ad53114e0137dc240c5f76f025..e9900b4f1c8e3cacb9735b3feb97cfaa6772b5a7 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -37,7 +37,7 @@ class gmshFace : public GFace {
   virtual SPoint2 geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t);
   virtual GPoint point(double par1, double par2) const; 
   virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const; 
-  virtual int containsPoint(const SPoint3 &pt) const;  
+  virtual bool containsPoint(const SPoint3 &pt) const;  
   virtual double getMetricEigenvalue(const SPoint2 &);  
   virtual SVector3 normal(const SPoint2 &param) const; 
   virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const; 
diff --git a/Geo/gmshSurface.cpp b/Geo/gmshSurface.cpp
index e236ae84696d150025d3dfe822cd51bf281a1407..357ef9a6adc7be26c636a5b97398bc7088b32c7a 100644
--- a/Geo/gmshSurface.cpp
+++ b/Geo/gmshSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshSurface.cpp,v 1.15 2008-05-04 08:31:14 geuzaine Exp $
+// $Id: gmshSurface.cpp,v 1.16 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -27,6 +27,30 @@
 
 std::map<int,gmshSurface*> gmshSurface::allGmshSurfaces;
 
+SPoint2 gmshSurface::parFromPoint(double x, double y, double z)
+{
+  Msg::Error("Parametric coordinate computation not implemented for this type of surface");
+  return SPoint2();
+}
+
+SVector3 gmshSurface::normal(const SPoint2 &param) const
+{
+  Msg::Error("Normal computation not implemented for this type of surface");
+  return SVector3();
+}
+
+Pair<SVector3, SVector3> gmshSurface::firstDer(const SPoint2 &param)
+{
+  Msg::Error("First derivative not implemented for this type of surface");
+  return Pair<SVector3, SVector3>();
+}
+
+double gmshSurface::getMetricEigenvalue(const SPoint2 &)
+{
+  Msg::Error("Metric eigenvalue not implemented for this type of surface");
+  return 0;
+}
+
 gmshSurface *gmshSphere::NewSphere(int iSphere, double x, double y, double z, double r)
 {
   gmshSphere *sph = new gmshSphere(x, y, z, r);
@@ -55,7 +79,6 @@ SPoint3 gmshSphere::point(double par1, double par2) const
   const double x = xc + r * sin(par2) * cos(par1);
   const double y = yc + r * sin(par2) * sin(par1);
   const double z = zc - r * cos(par2); 
-  // printf("%g %g - %g %g %g\n",par1,par2,x,y,z);
   return SPoint3(x, y, z);
 }
 
@@ -78,10 +101,10 @@ SPoint3 gmshPolarSphere::point(double parA, double parB) const
   //at the center of the sphere 
   //parA=2rx/(r+z) parB=2ry/(r+z)
   double rp2 = parA * parA + parB * parB;
-        SPoint3 p(2*parA/(1+rp2),2*parB/(1+rp2),(rp2-1)/(rp2+1));
-        p*=-r;
-        p+=o;
-        return p;
+  SPoint3 p(2*parA/(1+rp2),2*parB/(1+rp2),(rp2-1)/(rp2+1));
+  p *= -r;
+  p += o;
+  return p;
 }
 
 gmshSurface *gmshParametricSurface::NewParametricSurface(int iSurf, char *valX,
@@ -132,3 +155,9 @@ SPoint3 gmshParametricSurface::point(double par1, double par2) const
   return SPoint3(x, y, z);
 #endif
 }
+
+Range<double> gmshParametricSurface::parBounds(int i) const
+{
+  Msg::Error("Parameter bounds not available for parametric surface");
+  return Range<double>(0., 0.);
+}
diff --git a/Geo/gmshSurface.h b/Geo/gmshSurface.h
index 7bab38f55c435592b8bfb896a39ef7ad32e62adf..519710b00fee02cbbeb6a0ba7e7c14725497b2ca 100644
--- a/Geo/gmshSurface.h
+++ b/Geo/gmshSurface.h
@@ -59,12 +59,12 @@ public:
   virtual gmshSurface::gmshSurfaceType geomType() const = 0;
   virtual SPoint3 point(double par1, double par2) const = 0;
   virtual SPoint3 point(const SPoint2 &p) const { return point(p.x(), p.y()); }
-  virtual SPoint2 parFromPoint(double x, double y, double z) const = 0;
+  virtual SPoint2 parFromPoint(double x, double y, double z);
   // Return the normal to the face at the given parameter location.
-  virtual SVector3 normal(const SPoint2 &param) const = 0;
+  virtual SVector3 normal(const SPoint2 &param) const;
   // Return the first derivate of the face at the parameter location.
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const = 0;
-  virtual double getMetricEigenvalue(const SPoint2 &) { throw; }
+  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param);
+  virtual double getMetricEigenvalue(const SPoint2 &);
 };
 
 class gmshSphere : public gmshSurface
@@ -75,19 +75,13 @@ public:
   static gmshSurface *NewSphere(int _iSphere, double _x, double _y, double _z, double _r);
   virtual Range<double> parBounds(int i) const 
   { 
-    if(i == 0) return Range<double>(0., 2 * M_PI);
-    if(i == 1) return Range<double>(0., M_PI);
-    throw;
+    if(i == 0) 
+      return Range<double>(0., 2 * M_PI);
+    else
+      return Range<double>(0., M_PI);
   }
-  // Underlying geometric representation of this entity.
   virtual gmshSurface::gmshSurfaceType geomType() const { return gmshSurface::Sphere; }
   virtual SPoint3 point(double par1, double par2) const;
-  virtual SPoint2 parFromPoint(double x, double y, double z) const 
-  {
-    // 2 be done 
-    throw;
-  }
-  // Return the normal to the face at the given parameter location.
   virtual SVector3 normal(const SPoint2 &param) const
   {
     SPoint3  p1 = gmshSurface::point(param);
@@ -96,13 +90,8 @@ public:
     n.normalize();
     return n;
   }
-  // Return the first derivate of the face at the parameter location.
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const
-  {
-    // 2 be done
-    throw;
-  }  
 };
+
 class gmshPolarSphere : public gmshSurface
 {
   double r;
@@ -112,19 +101,13 @@ public:
   static gmshSurface *NewPolarSphere(int _iSphere, double _x, double _y, double _z, double _r);
   virtual Range<double> parBounds(int i) const 
   { 
-    if(i == 0) return Range<double>(-M_PI, M_PI);
-    if(i == 1) return Range<double>(-M_PI, M_PI);
-    throw;
+    if(i == 0)
+      return Range<double>(-M_PI, M_PI);
+    else
+      return Range<double>(-M_PI, M_PI);
   }
-  // Underlying geometric representation of this entity.
   virtual gmshSurface::gmshSurfaceType geomType() const { return gmshSurface::PolarSphere; }
   virtual SPoint3 point(double par1, double par2) const;
-  virtual SPoint2 parFromPoint(double x, double y, double z) const 
-  {
-    // 2 be done 
-    throw;
-  }
-  // Return the normal to the face at the given parameter location.
   virtual SVector3 normal(const SPoint2 &param) const
   {
     SPoint3  p1 = gmshSurface::point(param);
@@ -132,18 +115,11 @@ public:
     n.normalize();
     return n;
   }
-  // Return the first derivate of the face at the parameter location.
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const
-  {
-    // 2 be done
-    throw;
-  }  
   virtual double getMetricEigenvalue ( const SPoint2 &p)
   {
     double l = (4*r*r)/(4*r*r+p.x()*p.x()+p.y()*p.y());
     return l*l;
   }
-
 };
 
 class gmshParametricSurface : public gmshSurface
@@ -153,30 +129,12 @@ class gmshParametricSurface : public gmshSurface
   ~gmshParametricSurface();
 public:
   static gmshSurface *NewParametricSurface(int iSurf, char*, char*, char*);
-  virtual Range<double> parBounds(int i) const 
+  virtual Range<double> parBounds(int i) const;
+  virtual gmshSurface::gmshSurfaceType geomType() const 
   { 
-    throw;
+    return gmshSurface::ParametricSurface; 
   }
-  // Underlying geometric representation of this entity.
-  virtual gmshSurface::gmshSurfaceType geomType() const { return gmshSurface::ParametricSurface; }
   virtual SPoint3 point(double par1, double par2) const;
-  virtual SPoint2 parFromPoint(double x, double y, double z) const 
-  {
-    // 2 be done 
-    throw;
-  }
-  // Return the normal to the face at the given parameter location.
-  virtual SVector3 normal(const SPoint2 &param) const
-  {
-    throw;
-  }
-  // Return the first derivate of the face at the parameter location.
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const
-  {
-    // 2 be done
-    throw;
-  }  
 };
 
-
 #endif
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index c57d6d8e83d4b7f5646cb25cd3b0145b48d755e7..e34a09049460bbd819a29ffbbc3dc6abb93100ad 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.107 2008-06-10 08:37:34 remacle Exp $
+// $Id: BDS.cpp,v 1.108 2008-07-03 17:06:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -292,7 +292,10 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
   BDS_Point *p1 = find_point(num1);
   BDS_Point *p2 = find_point(num2);
   
-  if(!p1 || !p2) throw;
+  if(!p1 || !p2) {
+    Msg::Fatal("Could not find points %d or %d in BDS mesh", num1, num2);
+    return 0;
+  }
 
   Msg::Debug("edge %d %d has to be recovered", num1, num2);
   
@@ -428,8 +431,10 @@ BDS_Edge *BDS_Mesh::add_edge(int p1, int p2)
 
   BDS_Point *pp1 = find_point(p1);
   BDS_Point *pp2 = find_point(p2);
-  if(!pp1 || !pp2)
-    throw;
+  if(!pp1 || !pp2){
+    Msg::Fatal("Could not find points %d or %d in BDS mesh", p1, p2);
+    return 0;
+  }
   BDS_Edge *e = new BDS_Edge(pp1, pp2);
   edges.push_back(e);
   return e;
@@ -1230,8 +1235,10 @@ bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face
     pd[0] = u;
     pd[1] = v;
   }
-  else
-    throw;
+  else{
+    Msg::Error("Something wrong in move_point_parametric_quad");
+    return false;
+  }
   
   double ori_final1 = gmsh::orient2d(pa, pb, pc);
   double ori_final2 = gmsh::orient2d(pc, pd, pa);
@@ -1239,8 +1246,6 @@ bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face
   return ori_init1*ori_final1 > 0 && ori_init2*ori_final2 > 0 ;
 }
 
-
-
 // d^2_i = (x^2_i - x)^T M (x_i - x)  
 //       = M11 (x_i - x)^2 + 2 M21 (x_i-x)(y_i-y) + M22 (y_i-y)^2        
 
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 7b8193c66fa8a1b8743b7b2e5077a7e7c06c3798..5d41bebc91e8e98babf9936d92051de97a9c8718 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -25,16 +25,14 @@
 // points may know the normals to the surface they are classified on
 // default values are 0,0,1
 
-#include <stdio.h>
-#include <string>
 #include <set>
-#include <map>
 #include <vector>
 #include <algorithm>
 #include <list>
 #include <math.h>
 #include "GFace.h"
 #include "PView.h"
+#include "Message.h"
 
 class BDS_Edge;
 class BDS_Face;
@@ -56,11 +54,9 @@ void swap_config(BDS_Edge *e,
                  BDS_Point **p31, BDS_Point **p32, BDS_Point **p33,
                  BDS_Point **p41, BDS_Point **p42, BDS_Point **p43);
 
-
 class BDS_GeomEntity
 {
 public:
-
   int classif_tag;
   int classif_degree;
   inline bool operator < (const BDS_GeomEntity & other) const
@@ -70,21 +66,14 @@ public:
     if(classif_tag < other.classif_tag)return true;
     return false;
   }
-  BDS_GeomEntity(int a, int b)  
-    : classif_tag(a),classif_degree(b)
-  {
-  }
-  ~BDS_GeomEntity()  
-  {
-  }
+  BDS_GeomEntity(int a, int b) : classif_tag(a),classif_degree(b) {}
+  ~BDS_GeomEntity(){}
 };
 
-void print_face(BDS_Face *t);
-
 class BDS_Vector
 {
 public:
-  double x,y,z;
+  double x, y, z;
   bool operator < (const BDS_Vector &o) const
   {
     if(x - o.x  > t ) return true;
@@ -169,9 +158,9 @@ public:
 
 class BDS_Point 
 {
-  // the first size is the one dictated by the Background Mesh
-  // the second one is dictated by charecteristic lengths at points
-  // and is propagated
+  // the first size is the one dictated by the Background Mesh the
+  // second one is dictated by charecteristic lengths at points and is
+  // propagated
   double _lcBGM, _lcPTS;
 public:
   double X, Y, Z;
@@ -180,11 +169,9 @@ public:
   int iD;
   BDS_GeomEntity *g;
   std::list<BDS_Edge*> edges;
-
   // just a transition
   double &lcBGM() { return _lcBGM; }
   double &lc() { return _lcPTS; }
-  
   inline bool operator < (const BDS_Point & other) const
   {
     return iD < other.iD;
@@ -256,13 +243,14 @@ public:
   inline BDS_Face * otherFace(const BDS_Face *f) const
   {
     if(numfaces()!=2) {
-      printf("otherFace wrong, ony %d faces attached to edge %d %d\n",numfaces(),p1->iD,p2->iD);
-      throw;
+      Msg::Fatal("otherFace wrong, ony %d faces attached to edge %d %d",
+		 numfaces(), p1->iD, p2->iD);
+      return 0;
     }
     if(f == _faces[0]) return _faces[1];
     if(f == _faces[1]) return _faces[0];
-    printf("otherFace wrong : the edge does not belong to the face \n");
-    throw;
+    Msg::Fatal("otherFace wrong: the edge does not belong to the face");
+    return 0;
   }
   inline void del(BDS_Face *t)
   {
@@ -270,18 +258,14 @@ public:
                                 std::bind2nd(std::equal_to<BDS_Face*>(), t)), 
                  _faces.end());
   }
-  
   void oppositeof(BDS_Point * oface[2]) const; 
-  
   void update()
   {
     _length = sqrt((p1->X - p2->X) * (p1->X - p2->X) + 
                    (p1->Y - p2->Y) * (p1->Y - p2->Y) + 
                    (p1->Z - p2->Z) * (p1->Z - p2->Z));
   }
-
-  BDS_Edge(BDS_Point *A, BDS_Point *B)
-    : deleted(false), g(0)
+  BDS_Edge(BDS_Point *A, BDS_Point *B) : deleted(false), g(0)
   {         
     if(*A < *B){
       p1 = A;
@@ -304,28 +288,31 @@ public:
   BDS_Edge *e1, *e2, *e3, *e4;
   BDS_GeomEntity *g;
   inline int numEdges () const { return e4 ? 4 : 3; }
-  inline BDS_Edge *oppositeEdge (BDS_Point *p){
+  inline BDS_Edge *oppositeEdge (BDS_Point *p)
+  {
     if (e4){
-      printf("oppositeEdge to point %d cannot be applied to a quad\n",p->iD);
-      throw;
+      Msg::Fatal("oppositeEdge to point %d cannot be applied to a quad", p->iD);
+      return 0;
     }
-    if (e1->p1 != p && e1->p2 != p)return e1;
-    if (e2->p1 != p && e2->p2 != p)return e2;
-    if (e3->p1 != p && e3->p2 != p)return e3;
-    printf("point %d does not belong to this triangle\n",p->iD);
-    throw;
+    if (e1->p1 != p && e1->p2 != p) return e1;
+    if (e2->p1 != p && e2->p2 != p) return e2;
+    if (e3->p1 != p && e3->p2 != p) return e3;
+    Msg::Fatal("point %d does not belong to this triangle", p->iD);
+    return 0;
   }
-  inline BDS_Point *oppositeVertex (BDS_Edge *e){
+  inline BDS_Point *oppositeVertex (BDS_Edge *e)
+  {
     if (e4){
-      printf("oppositeVertex to edge %d %d cannot be applied to a quad\n",e->p1->iD,e->p2->iD);
-      throw;
+      Msg::Fatal("oppositeVertex to edge %d %d cannot be applied to a quad", 
+		 e->p1->iD, e->p2->iD);
+      return 0;
     }
-
-    if (e == e1)return e2->commonvertex(e3);
-    if (e == e2)return e1->commonvertex(e3);
-    if (e == e3)return e1->commonvertex(e2);
-    printf("edge  %d %d does not belong to this triangle\n",e->p1->iD,e->p2->iD);
-    throw;
+    if (e == e1) return e2->commonvertex(e3);
+    if (e == e2) return e1->commonvertex(e3);
+    if (e == e3) return e1->commonvertex(e2);
+    Msg::Fatal("edge  %d %d does not belong to this triangle",
+	       e->p1->iD, e->p2->iD);
+    return 0;
   }
   inline void getNodes(BDS_Point *n[4]) const
   {
@@ -342,7 +329,6 @@ public:
       n[3] = e3->commonvertex(e4);
     }
   }
-  
   BDS_Face(BDS_Edge *A, BDS_Edge *B, BDS_Edge *C,BDS_Edge *D = 0)
     : deleted(false), e1(A), e2(B), e3(C), e4(D), g(0)
   {     
@@ -506,7 +492,6 @@ public:
 
 void outputScalarField(std::list<BDS_Face*> t, const char *fn, int param, GFace *gf=0);
 void recur_tag(BDS_Face *t, BDS_GeomEntity *g);
-
 int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
                        double x3, double y3, double x4, double y4,
 		       double x[2]);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index ec7e51b86c5c0cd1c492230eea5b2efae6f1ecf0..20dad12205cc54d9f6df85d20a908265137130f7 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.145 2008-06-05 11:52:49 samtech Exp $
+// $Id: Generator.cpp,v 1.146 2008-07-03 17:06:03 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -413,7 +413,7 @@ void GenerateMesh(GModel *m, int ask)
   // Orient the surface mesh so that it matches the geometry
   if(m->getMeshStatus() >= 2)
     std::for_each(m->firstFace(), m->lastFace(), orientMeshGFace());
-  
+
   // Optimize quality of 3D tet mesh
   if(m->getMeshStatus() == 3){
     for(int i = 0; i < std::max(CTX.mesh.optimize, CTX.mesh.optimize_netgen); i++){
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 5d5096205bee412c4f6aefb4f288d9cb74c711d3..a67d10c4d16fe41018bf4da3d4b552904c8932ff 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.31 2008-06-20 12:15:44 remacle Exp $
+// $Id: HighOrder.cpp,v 1.32 2008-07-03 17:06:03 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -503,7 +503,6 @@ void getFaceVertices(GFace *gf,
           }
         }
       }
-      else throw; // not tri or quad
     }
   }  
 }
@@ -603,7 +602,6 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
           }
         }
       }
-      else throw; // not tri or quad
     }
   }
 }
@@ -646,7 +644,6 @@ void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
           }
         }
       }
-      else throw; // not tri or quad
     }
   }
 }
@@ -1063,8 +1060,14 @@ void optimizeNodeLocations(GFace *gf, smoothVertexDataHON &vdN, double eps = .2)
   if(!vdN.v.size()) return;
   double uv[20];
   for (unsigned int i = 0; i < vdN.v.size(); i++){
-    if (!vdN.v[i]->getParameter(0, uv[2 * i])) throw;
-    if (!vdN.v[i]->getParameter(1, uv[2 * i + 1])) throw;
+    if (!vdN.v[i]->getParameter(0, uv[2 * i])){
+      Msg::Error("Node location optimization failed");
+      return;
+    }
+    if (!vdN.v[i]->getParameter(1, uv[2 * i + 1])){
+      Msg::Error("Node location optimization failed");
+      return;
+    }
   }
 
   double F = -smooth_obj_HighOrderN(uv, &vdN);
@@ -1420,6 +1423,11 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
   //   edges (creating 8-node quads, 20-node hexas, etc., instead of
   //   9-node quads, 27-node hexas, etc.)
 
+#if !defined(HAVE_GSL)
+  Msg::Error("High order mesh generation requires the GSL");
+  return;
+#endif
+
   int nPts = order - 1;
 
   Msg::StatusBar(1, true, "Meshing second order...");
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 096bb0bf07cd291b29e0fa9946dd6b1dfbbbedbe..b3ced4eb27ef8591fea6c6c8901302e0ea106e81 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.226 2008-06-20 05:51:37 geuzaine Exp $
+# $Id: Makefile,v 1.227 2008-07-03 17:06:03 geuzaine Exp $
 #
 # Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 #
@@ -215,12 +215,12 @@ meshGFaceDelaunayInsertion.o: meshGFaceDelaunayInsertion.cpp BDS.h \
   ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
-  ../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h BackgroundMesh.h \
-  meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
-  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
-  meshGFaceOptimize.h meshGFace.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Common/Message.h ../Common/Context.h
+  ../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/Message.h \
+  BackgroundMesh.h meshGFaceDelaunayInsertion.h ../Geo/MElement.h \
+  ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
+  ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
+  ../Geo/MVertex.h ../Geo/SVector3.h meshGFaceOptimize.h meshGFace.h \
+  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Common/Context.h
 meshGFaceOptimize.o: meshGFaceOptimize.cpp meshGFaceOptimize.h \
   ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
   ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
@@ -231,7 +231,7 @@ meshGFaceOptimize.o: meshGFaceOptimize.cpp meshGFaceOptimize.h \
   ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
-  ../Geo/Pair.h BackgroundMesh.h
+  ../Geo/Pair.h BackgroundMesh.h ../Common/Message.h
 meshGFaceQuadrilateralize.o: meshGFaceQuadrilateralize.cpp \
   meshGFaceQuadrilateralize.h ../Common/Message.h ../Numeric/Numeric.h \
   ../Numeric/NumericEmbedded.h meshGFaceDelaunayInsertion.h \
@@ -346,11 +346,11 @@ qualityMeasures.o: qualityMeasures.cpp qualityMeasures.h BDS.h \
   ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
-  ../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Geo/MVertex.h \
-  ../Geo/SPoint3.h ../Geo/MElement.h ../Common/GmshDefines.h \
-  ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
-  ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h
+  ../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/Message.h \
+  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MElement.h \
+  ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/MEdge.h \
+  ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+  ../Geo/SVector3.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h
 BoundaryLayers.o: BoundaryLayers.cpp ../Geo/GModel.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index a0d7c9d0b447cb73a82bb6ace7a721da727405c8..359523254c6ddc27946354f51b70ddd100f3d364 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.138 2008-07-01 14:24:07 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.139 2008-07-03 17:06:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -196,9 +196,15 @@ void computeEdgeLoops(const GFace *gf,
       v_start = start;
     }
     else{
-      if(it == edges.end ()) throw;
+      if(it == edges.end ()){
+	Msg::Error("Something wrong in edge loop computation");
+	return;
+      }
       v_start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
-      if(v_start != v_end) throw;
+      if(v_start != v_end){
+	Msg::Error("Something wrong in edge loop computation");
+	return;
+      }
       v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
     }
     all_mvertices.push_back(v_start->mesh_vertices[0]);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index ab3db5a635fd713dcdc5f48720d393142c030eaf..b0bbd1f784187bad54abe695db314801fbb7beda 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.27 2008-06-03 12:43:42 remacle Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.28 2008-07-03 17:06:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -557,7 +557,10 @@ static void insertAPoint(GFace *gf,
 			 MTri3 *worst = 0){
   if (worst){
     it = AllTris.find(worst);
-    if (worst != *it)throw;
+    if (worst != *it){
+      Msg::Error("Could not insert point");
+      return;
+    }
   }
   else worst = *it;
 
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 345efb6fb3cb97de6fef23ebc2865512b4649464..453d52305e8eb32a89c39e2a235163b3145dd798 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceOptimize.cpp,v 1.15 2008-03-26 20:32:40 remacle Exp $
+// $Id: meshGFaceOptimize.cpp,v 1.16 2008-07-03 17:06:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -27,6 +27,7 @@
 #include "MVertex.h"
 #include "MElement.h"
 #include "BackgroundMesh.h"
+#include "Message.h"
 
 static void setLcsInit(MTriangle *t, std::map<MVertex*, double> &vSizes)
 {
@@ -325,7 +326,8 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
     }
     break;
   default :
-    throw;
+    Msg::Error("Unknown swapping criterion");
+    return false;
   }
 
   std::list<MTri3*> cavity;
@@ -453,7 +455,8 @@ bool gmshEdgeSplit(const double lMax, MTri3 *t1, GFace *gf, int iLocalEdge,
   case SPCR_ALLWAYS:
     break;
   default :
-    throw;
+    Msg::Error("Unknown swapping criterion");
+    return false;
   }
 
   gf->mesh_vertices.push_back(vnew);
@@ -584,7 +587,10 @@ bool gmshBuildVertexCavity(MTri3 *t, int iLocalVertex, MVertex **v1,
           return true;
         }
         if (!t) return false;
-        if (t->isDeleted()){ printf("weird\n"); throw; }  
+        if (t->isDeleted()){ 
+	  Msg::Error("Impossible to build vertex cavity");
+	  return false;
+	}  
         cavity.push_back(t);
         for (int j = 0; j < 3; j++){
           if (t->tri()->getVertex(j) !=lastinring && t->tri()->getVertex(j) != *v1){
@@ -596,7 +602,10 @@ bool gmshBuildVertexCavity(MTri3 *t, int iLocalVertex, MVertex **v1,
         break;
       }
     }
-    if (iEdge == -1) { printf("not found\n"); throw; }
+    if (iEdge == -1) {
+      Msg::Error("Impossible to build vertex cavity");
+      return false;
+    }
   }
 }
 
diff --git a/Mesh/meshGFaceQuadrilateralize.cpp b/Mesh/meshGFaceQuadrilateralize.cpp
index ac1df83780efcfe0ca981298b11c3bb13cc17146..fef2f65da3a214df4867c88c38f880119045ca9e 100644
--- a/Mesh/meshGFaceQuadrilateralize.cpp
+++ b/Mesh/meshGFaceQuadrilateralize.cpp
@@ -96,7 +96,7 @@ void gmshEdgeFront::updateStatus (BDS_Edge *e){
       return;
     }
   }
-  throw;
+  Msg::Error("Something wrong in updateStatus");
 }
 
 SVector3 norm_edge (BDS_Point *p1, BDS_Point *p2){
@@ -172,8 +172,10 @@ SVector3 gmshEdgeFront::normal (BDS_Edge*e) const{
     p3 = t->e1->commonvertex(t->e3);
   else if (e == t->e3)
     p3 = t->e2->commonvertex(t->e1);
-  else 
-    throw;
+  else {
+    Msg::Error("Could not compute fron normal");
+    return SVector3();
+  }
   
   SVector3 t1 (p2->X-p1->X,p2->Y-p1->Y,p2->Z-p1->Z);
   SVector3 t2 (p3->X-p1->X,p3->Y-p1->Y,p3->Z-p1->Z);
@@ -209,9 +211,7 @@ void gmshEdgeFront::getFrontEdges (BDS_Point *p, eiter & it1, eiter & it2) const
       if ( it2 != edges.end()) return;
     }
   }  
-  printf("point %d is in the front but has only %d edges\n",p->iD,count);
-
-  throw;
+  Msg::Error("point %d is in the front but has only %d edges\n",p->iD,count);
 }
 
 void gmshEdgeFront::initiate () {
@@ -513,7 +513,8 @@ bool gmshEdgeFront::emptyFront (int tag){
     //    printf("strange\n");
     break;
   default:
-    throw;
+    Msg::Error("Unknown case in emptyFront");
+    return false;
   }
 
   if ( fe1.size() || fe2.size() || !left || !right || !formQuad (e, left, right)){
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index 431eb6c3c3be2c1e0cb55d10a95ecfff4af95ae5..d038e9a55f8036d7c9dbc588377f46fc1ec07f90 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionLocalMeshMod.cpp,v 1.13 2008-05-04 08:31:16 geuzaine Exp $
+// $Id: meshGRegionLocalMeshMod.cpp,v 1.14 2008-07-03 17:06:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -99,10 +99,10 @@ bool gmshBuildEdgeCavity(MTet4 *t,
     else if (faces[iFace2][0] == edges[5-iLocalEdge][K] ||
              faces[iFace2][1] == edges[5-iLocalEdge][K] ||
              faces[iFace2][2] == edges[5-iLocalEdge][K] ) iFace = iFace2;
-    else { Msg::Error("Error of connexion"); throw; }
+    else { Msg::Error("Error of connexion"); return false; }
     t=t->getNeigh(iFace);
     if (!t) return false;
-    if (t->isDeleted()) throw;
+    if (t->isDeleted()){ Msg::Error("Weird!!"); return false; }
     if (t == cavity[0]) break;
     ring.push_back(lastinring);    
     cavity.push_back(t);
@@ -117,7 +117,7 @@ bool gmshBuildEdgeCavity(MTet4 *t,
     }  
     if (iLocalEdge == -1){
       Msg::Error("loc = %d", iLocalEdge);
-      throw;
+      return false;
     }
   }
   computeNeighboringTetsOfACavity (cavity,outside);
@@ -384,7 +384,10 @@ bool gmshFaceSwap(std::vector<MTet4 *> &newTets,
       break;
     }
   }
-  if (!v2) throw;
+  if (!v2){
+    Msg::Error("Impossible to swap face");
+    return false;
+  }
 
   // printf("%p %p -- %p %p %p\n",v1,v2,f1,f2,f3);
 
@@ -573,7 +576,10 @@ bool gmshCollapseVertex(std::vector<MTet4 *> &newTets,
                         const gmshLocalMeshModAction action,
                         double *minQual)
 {
-  if(t->isDeleted()) throw;
+  if(t->isDeleted()){
+    Msg::Error("Impossible to collapse vertex");
+    return false;
+  }
 
   MVertex *v = t->tet()->getVertex(iVertex);
   MVertex *tg = t->tet()->getVertex(iTarget);
@@ -613,7 +619,10 @@ bool gmshCollapseVertex(std::vector<MTet4 *> &newTets,
   
   double worstAfter = 1.0;
   double newQuals[2000];
-  if (toUpdate.size() >= 2000) throw;
+  if (toUpdate.size() >= 2000){
+    Msg::Error("Impossible to collapse vertex");
+    return false;
+  }
   for (unsigned int i = 0; i < toUpdate.size(); i++){
     double vv;
     newQuals[i] = qmTet(toUpdate[i]->tet(),cr,&vv);
@@ -659,7 +668,10 @@ bool gmshSmoothVertex(MTet4 *t,
                       int iVertex,
                       const gmshQualityMeasure4Tet &cr)
 {
-  if(t->isDeleted()) throw;
+  if(t->isDeleted()){
+    Msg::Error("Impossible to collapse vertex");
+    return false;
+  }
   if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3) return false;
 
   std::vector<MTet4*> cavity;
@@ -702,7 +714,10 @@ bool gmshSmoothVertex(MTet4 *t,
   t->tet()->getVertex(iVertex)->z() = zcg;
   double worstAfter = 1.0;
   double newQuals[2000];
-  if (cavity.size() >= 2000) throw;
+  if (cavity.size() >= 2000){
+    Msg::Error("Impossible to smooth vertex");
+    return false;
+  }
   for (unsigned int i = 0; i < cavity.size(); i++){
     double volume;
     newQuals[i] = qmTet(cavity[i]->tet(),cr,&volume);
@@ -816,7 +831,10 @@ bool gmshSmoothVertexOptimize(MTet4 *t,
   t->tet()->getVertex(iVertex)->z() = zopti;
 
   double newQuals[2000];
-  if(vd.ts.size() >= 2000) throw;
+  if(vd.ts.size() >= 2000){
+    Msg::Error("Impossible to smooth vertex");
+    return false;
+  }
   for(unsigned int i = 0; i < vd.ts.size(); i++){
     double volume;
     newQuals[i] = qmTet(vd.ts[i]->tet(), cr, &volume);
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index d77fe61f93e0152ce11dd0818131dab35d9fbd09..b0fb0e450d960d0e8c802d54fc9b839152a9e7f1 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -1,4 +1,4 @@
-// $Id: qualityMeasures.cpp,v 1.13 2008-03-20 11:44:09 geuzaine Exp $
+// $Id: qualityMeasures.cpp,v 1.14 2008-07-03 17:06:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -24,6 +24,7 @@
 #include "MVertex.h"
 #include "MElement.h"
 #include "Numeric.h"
+#include "Message.h"
 
 double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3, 
                   const gmshQualityMeasure4Triangle &cr)
@@ -90,7 +91,8 @@ double qmTriangle(const double &xa, const double &ya, const double &za,
     }
     break;
   default:
-    throw;
+    Msg::Error("Unknown quality measure");
+    return 0.;
   }  
 
   return quality;
@@ -181,6 +183,7 @@ double qmTet(const double &x1, const double &y1, const double &z1,
     }
     break;
   default:
-    throw;
+    Msg::Error("Unknown quality measure");
+    return 0.;
   }
 }
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index 7cea8abaa66149a494e6d8630e0ee2c949fe3d8d..bebaa356a197e654f3067912875f663d8182f0b8 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -1,4 +1,4 @@
-// $Id: FunctionSpace.cpp,v 1.7 2008-06-27 08:10:07 koen Exp $
+// $Id: FunctionSpace.cpp,v 1.8 2008-07-03 17:06:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -21,6 +21,7 @@
 
 #include "FunctionSpace.h"
 #include "GmshDefines.h"
+#include "Message.h"
 
 Double_Matrix generatePascalTriangle(int order)
 {
@@ -503,8 +504,10 @@ Double_Matrix gmshGeneratePointsTriangle(int order, bool serendip)
 Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial,
                                                    const Double_Matrix& point) 
 {
-  if (monomial.size1() != point.size1()) throw;
-  if (monomial.size2() != point.size2()) throw;
+  if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
+    Msg::Fatal("Wrong sizes for Lagrange coefficients generation");
+    return Double_Matrix(1, 1);
+  }
   
   int ndofs = monomial.size1();
   int dim   = monomial.size2();
@@ -523,7 +526,10 @@ Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial
   
   double det = Vandermonde.determinant();
 
-  if (det == 0.0) throw;
+  if (det == 0.0){
+    Msg::Fatal("Vandermonde matrix has zero determinant!?");
+    return Double_Matrix(1, 1);
+  }
 
   Double_Matrix coefficient(ndofs, ndofs);
   
@@ -612,7 +618,10 @@ const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(5, false);
     break;
   default :
-    throw;
+    Msg::Error("Unknown function space %d: reverting to TET_4", tag);
+    F.monomials = generatePascalTetrahedron(1);
+    F.points =    gmshGeneratePointsTetrahedron(1, false);
+    break;
   }  
   F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
   fs.insert(std::make_pair(tag, F));
diff --git a/Numeric/GaussLegendre1D.h b/Numeric/GaussLegendre1D.h
index a209e8a8fd579f4cca4665f80ed7e39272811e36..d1124cee82ae399664462d9a4c9ae5e1ee279d44 100644
--- a/Numeric/GaussLegendre1D.h
+++ b/Numeric/GaussLegendre1D.h
@@ -214,7 +214,9 @@ inline void gmshGaussLegendre1D(int nbQuadPoints, double **t, double **w)
     *w = _GL_wt16;
     break;
   default :
-    throw;
+    *t = 0;
+    *w = 0;
+    break;
   }
 }
 
diff --git a/Numeric/Makefile b/Numeric/Makefile
index 86f342170cc48c2feb809ca048be5929e0e5d902..6faed686d00b2d8c4f69d8ef2c4bce59a754f60c 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.56 2008-06-07 17:20:48 geuzaine Exp $
+# $Id: Makefile,v 1.57 2008-07-03 17:06:04 geuzaine Exp $
 #
 # Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 #
@@ -67,7 +67,7 @@ NumericEmbedded.o: NumericEmbedded.cpp NumericEmbedded.h \
   ../Common/Message.h
 EigSolve.o: EigSolve.cpp
 FunctionSpace.o: FunctionSpace.cpp FunctionSpace.h ../Common/GmshMatrix.h \
-  ../Common/GmshDefines.h
+  ../Common/GmshDefines.h ../Common/Message.h
 gmsh_predicates.o: gmsh_predicates.cpp
 gsl_newt.o: gsl_newt.cpp ../Common/Message.h Numeric.h NumericEmbedded.h
 gsl_min.o: gsl_min.cpp ../Common/Message.h Numeric.h NumericEmbedded.h
diff --git a/utils/misc/driverOCC.cpp b/utils/misc/driverOCC.cpp
index 22cd22bd29baff53a5c8935264d9ab450db3d57b..3cc2f1b440e65b0fc899f6463f59a338be5ea13d 100644
--- a/utils/misc/driverOCC.cpp
+++ b/utils/misc/driverOCC.cpp
@@ -13,6 +13,15 @@
 #include <gmsh/GModel.h>
 #include <gmsh/MElement.h>
 
+class mymsg : public GmshMessage{
+public:
+  void operator()(std::string level, std::string msg)
+  {
+    printf("level=%s msg=%s\n", level.c_str(), msg.c_str());
+    if(level == "Fatal") throw "Fatal error in Gmsh";
+  }
+};
+
 int main(int argc, char **argv)
 {
   // create an OCC shape (by loading it from a brep file)
@@ -21,12 +30,20 @@ int main(int argc, char **argv)
   BRepTools::Read(shape, argv[1], builder);
   BRepTools::Clean(shape);
 
-  // import the shape in gmsh and mesh it
+  // initialize gmsh and set a message callback
   GmshInitialize(argc, argv);
+  mymsg c;
+  GmshSetMessageHandler(&c);
 
+  // create a model, import the shape, and mesh it
   GModel m;
   m.importOCCShape((void*)&shape, 0);
-  m.mesh(2);
+  try{
+    m.mesh(2);
+  }
+  catch(...){
+    printf("Unrecoverable error in gmsh--aborting mesh!\n");
+  }
 
   for(GModel::fiter it = m.firstFace(); it != m.lastFace(); ++it){
     GFace *f = *it;