diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 8ab9c5ce025eb5a0675b1fe67f2fa5f65e7b258e..5ecf14f61555d78744684b3279918395bfd05c90 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdge.cpp,v 1.34 2008-01-21 22:16:04 geuzaine Exp $
+// $Id: GEdge.cpp,v 1.35 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -26,6 +26,12 @@
 #include "MElement.h"
 #include "GmshDefines.h"
 
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
+#include "Message.h"
+#endif
+
 GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
   : GEntity(model, tag), v0(_v0), v1(_v1)
 {
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index d6401b31a59753e91baea409420d75d206deeaa4..766b2a1ae1516dc5521eaece6bce46729926d3af 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdgeLoop.cpp,v 1.8 2008-01-20 11:19:27 geuzaine Exp $
+// $Id: GEdgeLoop.cpp,v 1.9 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -21,7 +21,12 @@
 
 #include <algorithm>
 #include "GEdgeLoop.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
 #include "Message.h"
+#endif
 
 void GEdgeSigned::print() const
 {
diff --git a/Geo/GEntity.cpp b/Geo/GEntity.cpp
index ca76a0ac45f3f9efd1cf6baa67f8ecec0458906c..654f0d0b9cc41ca8e8a97ce96afcdced152c8930 100644
--- a/Geo/GEntity.cpp
+++ b/Geo/GEntity.cpp
@@ -1,4 +1,4 @@
-// $Id: GEntity.cpp,v 1.15 2007-09-22 20:35:18 geuzaine Exp $
+// $Id: GEntity.cpp,v 1.16 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -20,8 +20,13 @@
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GEntity.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
 #include "VertexArray.h"
 #include "Context.h"
+#endif
 
 extern Context_T CTX;
 
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 67600f5ab2f2bee6eb312e773639547700ab9304..ec5ee8fc70f0b4c74c50ef739fb9428639be324d 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.41 2008-01-21 22:16:04 geuzaine Exp $
+// $Id: GFace.cpp,v 1.42 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -23,17 +23,21 @@
 #include "GFace.h"
 #include "GEdge.h"
 #include "MElement.h"
-#include "Message.h"
-#include "Numeric.h"
-#include "Context.h"
 
-#if defined(HAVE_GSL)
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
+#if defined(HAVE_GMSH_EMBEDDED)
+#  include "GmshEmbedded.h"
 #else
-#define NRANSI
-#include "nrutil.h"
+#  include "Message.h"
+#  include "Numeric.h"
+#  include "Context.h"
+#  if defined(HAVE_GSL)
+#    include <gsl/gsl_vector.h>
+#    include <gsl/gsl_linalg.h>
+#  else
+#    define NRANSI
+#    include "nrutil.h"
 void dsvdcmp(double **a, int m, int n, double w[], double **v);
+#  endif
 #endif
 
 extern Context_T CTX;
@@ -194,6 +198,7 @@ void GFace::computeMeanPlane(const std::vector<MVertex*> &points)
 
 void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
 {
+#if !defined(HAVE_GMSH_EMBEDDED)
   // The concept of a mean plane computed in the sense of least
   // squares is fine for plane surfaces(!), but not really the best
   // one for non-plane surfaces. Indeed, imagine a quarter of a circle
@@ -380,6 +385,7 @@ end:
       }
     }
   }
+#endif
 }
 
 void GFace::getMeanPlaneData(double VX[3], double VY[3], 
@@ -443,6 +449,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
 		    double &U, double &V, const double relax,
 		    const bool onSurface) const
 {
+#if !defined(HAVE_GMSH_EMBEDDED)
   const double Precision = 1.e-8;
   const int MaxIter = 25;
   const int NumInitGuess = 11;
@@ -519,6 +526,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
     Msg(INFO, "point %g %g %g : Relaxation factor = %g",X,Y,Z, 0.75 * relax);
     XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
   }  
+#endif
 }
 
 SPoint2 GFace::parFromPoint(const SPoint3 &p) const
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 0ef09179a51599d882ba347df1811bea57999c0e..abc4cd784a8f8dd671fd0bb2bd8835ba23ea7bfe 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.54 2008-01-19 23:04:12 geuzaine Exp $
+// $Id: GModel.cpp,v 1.55 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -20,12 +20,17 @@
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GModel.h"
-#include "gmshSurface.h"
 #include "MElement.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
+#include "Message.h"
+#include "gmshSurface.h"
 #include "Field.h"
 #include "BackgroundMesh.h"
-#include "Message.h"
 #include "Context.h"
+#endif
 
 extern Context_T CTX;
 
@@ -89,9 +94,11 @@ void GModel::destroy()
 
   MVertex::resetGlobalNumber();
   MElement::resetGlobalNumber();
+#if !defined(HAVE_GMSH_EMBEDDED)
   gmshSurface::reset();
   fields.reset();
   BGMReset();
+#endif
 }
 
 GRegion *GModel::regionByTag(int n) const
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 416e16b3e99aa549ac2e6ec54ac451baf157f648..ef7a5cd57006fcf7df5d7109a85c4041bde25e5a 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.27 2008-01-21 23:28:52 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.28 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -23,15 +23,20 @@
 #include <cstring>
 #include <string>
 
-#include "Message.h"
-#include "GmshDefines.h"
 #include "GModel.h"
+#include "GmshDefines.h"
+#include "MElement.h"
+#include "SBoundingBox3d.h"
 #include "discreteRegion.h"
 #include "discreteFace.h"
 #include "discreteEdge.h"
 #include "discreteVertex.h"
-#include "MElement.h"
-#include "SBoundingBox3d.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
+#include "Message.h"
+#endif
 
 static void swapBytes(char *array, int size, int n)
 {
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 8a429ef30a4c8cd021b9d1efc53cdba2f949cbda..c2ff22677779367ffe5a4b94f56489a005355bb7 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: GRegion.cpp,v 1.19 2008-01-19 22:06:01 geuzaine Exp $
+// $Id: GRegion.cpp,v 1.20 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -24,6 +24,12 @@
 #include "GFace.h"
 #include "MElement.h"
 
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
+#include "Message.h"
+#endif
+
 GRegion::GRegion(GModel *model, int tag) : GEntity (model, tag)
 {
   resetMeshAttributes();
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index a2480a6b37af3a27c99efc8cc408a47b8bdb4c38..06342e3dd9efe6741e313a66c9c41f18f1905da0 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: GVertex.cpp,v 1.15 2008-01-20 10:10:41 geuzaine Exp $
+// $Id: GVertex.cpp,v 1.16 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -24,7 +24,12 @@
 #include "GVertex.h"
 #include "GFace.h"
 #include "MVertex.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
 #include "Message.h"
+#endif
 
 GVertex::GVertex(GModel *m, int tag, double ms) : GEntity (m, tag), meshSize (ms) 
 {
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index c0f55b7ac6bad60dc8820cb32727974f81814a72..6a4c35c79ecbe7b9b331f6e0887bdad007f0d314 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.49 2008-01-22 13:48:48 geuzaine Exp $
+// $Id: MElement.cpp,v 1.50 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -23,10 +23,15 @@
 #include "MElement.h"
 #include "GEntity.h"
 #include "GFace.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
 #include "Numeric.h"
 #include "Message.h"
 #include "Context.h"
 #include "qualityMeasures.h"
+#endif
 
 extern Context_T CTX;
 
@@ -99,19 +104,31 @@ double MElement::rhoShapeMeasure()
 
 double MTriangle::gammaShapeMeasure()
 {
+#if defined(HAVE_GMSH_EMBEDDED)
+  return 0.;
+#else
   return qmTriangle(this, QMTRI_RHO);
+#endif
 }
 
 double MTetrahedron::gammaShapeMeasure()
 {
+#if defined(HAVE_GMSH_EMBEDDED)
+  return 0.;
+#else
   double vol;
   return qmTet(this, QMTET_2, &vol);
+#endif
 }
 
 double MTetrahedron::etaShapeMeasure()
 {
+#if defined(HAVE_GMSH_EMBEDDED)
+  return 0.;
+#else
   double vol;
   return qmTet(this, QMTET_3, &vol);
+#endif
 }
 
 void MTetrahedron::getMat(double mat[3][3])
diff --git a/Geo/MFace.cpp b/Geo/MFace.cpp
index 47a87f7810a429ca3ebee4991c6a22890fa764f0..f3b0de99b0f960b9705ee6f420c70de40aae7e35 100644
--- a/Geo/MFace.cpp
+++ b/Geo/MFace.cpp
@@ -18,8 +18,13 @@
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "MFace.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
 #include "Numeric.h"
 #include "Context.h"
+#endif
 
 extern Context_T CTX;
 
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 0c7532683e6739bae793cdbd913d0ab6fee2440d..a8582a09cab1f8ba1adf78315c2652ae9455b27e 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: MVertex.cpp,v 1.17 2007-12-03 15:17:40 remacle Exp $
+// $Id: MVertex.cpp,v 1.18 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -25,6 +25,12 @@
 #include "GEdge.h"
 #include "GFace.h"
 
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
+#include "Message.h"
+#endif
+
 int MVertex::_globalNum = 0;
 double MVertexLessThanLexicographic::tolerance = 1.e-6;
 
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 64d9fcb8d3f79df54764e92dde7752a5cf4ffbf4..66e8b5997eac2c3b2b4189472288e3615af3bfc9 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -4,7 +4,7 @@
 #include "GModel.h"
 #include "GEdge.h"
 
-#if !defined(HAVE_NO_GEO)
+#if !defined(HAVE_GMSH_EMBEDDED)
 #include "Geo.h"
 #endif
 
@@ -12,7 +12,7 @@ class discreteEdge : public GEdge {
  public:
   discreteEdge(GModel *model, int num) : GEdge(model, num, 0, 0) 
   {
-#if !defined(HAVE_NO_GEO)
+#if !defined(HAVE_GMSH_EMBEDDED)
     Curve *c = Create_Curve(num, MSH_SEGM_DISCRETE, 0, 0, 0, -1, -1, 0., 1.);
     Tree_Add(model->getGEOInternals()->Curves, &c);
     CreateReversedCurve(c);
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index cca31e4c63ae555006139ceb22e867d5585cb6b2..af9a3c40ab80465467b70fb0aba59639d22a6117 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -4,7 +4,7 @@
 #include "GModel.h"
 #include "GFace.h"
 
-#if !defined(HAVE_NO_GEO)
+#if !defined(HAVE_GMSH_EMBEDDED)
 #include "Geo.h"
 #endif
 
@@ -12,7 +12,7 @@ class discreteFace : public GFace {
  public:
   discreteFace(GModel *model, int num) : GFace(model, num)
   {
-#if !defined(HAVE_NO_GEO)
+#if !defined(HAVE_GMSH_EMBEDDED)
     Surface *s = Create_Surface(num, MSH_SURF_DISCRETE);
     Tree_Add(model->getGEOInternals()->Surfaces, &s);
 #endif
diff --git a/Geo/discreteRegion.h b/Geo/discreteRegion.h
index 71120f6f72d6478cd99ff06ff4cfa2ce79bd4f92..af2a58485ac1c7830383881e35f44ebada361da5 100644
--- a/Geo/discreteRegion.h
+++ b/Geo/discreteRegion.h
@@ -3,7 +3,7 @@
 
 #include "GRegion.h"
 
-#if !defined(HAVE_NO_GEO)
+#if !defined(HAVE_GMSH_EMBEDDED)
 #include "Geo.h"
 #endif
 
@@ -11,7 +11,7 @@ class discreteRegion : public GRegion {
  public:
   discreteRegion(GModel *model, int num) : GRegion(model, num)
   {
-#if !defined(HAVE_NO_GEO)
+#if !defined(HAVE_GMSH_EMBEDDED)
     ::Volume *v = Create_Volume(num, MSH_VOLUME_DISCRETE);
     Tree_Add(model->getGEOInternals()->Volumes, &v);
 #endif
diff --git a/Makefile b/Makefile
index cf0ad2906037ce2fc0ddc3a4fe643e77448cf928..f9f9d548c94e7e0eab6d7b1670ceb80b652b4686 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.458 2008-01-21 20:42:16 geuzaine Exp $
+# $Id: Makefile,v 1.459 2008-01-22 16:47:10 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -33,6 +33,14 @@ GMSH_SHORT_LICENSE = "GNU General Public License"
 GMSH_VERSION_FILE = Common/GmshVersion.h
 GMSH_DATE = `date "+%Y%m%d"`
 
+GMSH_API = Geo/GModel.h\
+           Geo/GEntity.h Geo/GPoint.h\
+           Geo/GVertex.h Geo/GEdge.h Geo/GEdgeLoop.h Geo/GFace.h Geo/GRegion.h\
+           Geo/MVertex.h Geo/MEdge.h Geo/MFace.h Geo/MElement.h\
+           Geo/SPoint2.h Geo/SPoint3.h Geo/SVector3.h Geo/SBoundingBox3d.h\
+           Geo/Pair.h Geo/Range.h\
+           Common/GmshDefines.h Common/GmshVersion.h
+
 all: link
 
 link: compile
@@ -66,16 +74,17 @@ lib: compile
 install-lib: lib
 	mkdir -p ${includedir}/gmsh
 	rm -f ${includedir}/gmsh/*
-	cp -f Geo/GModel.h\
-                Geo/GEntity.h Geo/GPoint.h\
-                Geo/GVertex.h Geo/GEdge.h Geo/GEdgeLoop.h Geo/GFace.h Geo/GRegion.h\
-              Geo/MVertex.h Geo/MEdge.h Geo/MFace.h Geo/MElement.h\
-              Geo/SPoint2.h Geo/SPoint3.h Geo/SVector3.h Geo/SBoundingBox3d.h\
-              Geo/Pair.h Geo/Range.h\
-              Common/GmshDefines.h Common/GmshVersion.h\
-          ${includedir}/gmsh
+	cp -f ${GMSH_API} ${includedir}/gmsh
 	cp -f bin/libGmsh${LIBEXT} ${libdir}
 
+embed:
+	cp -f ${GMSH_API} Geo/discrete*.h Numeric/NumericEmbedded.h utils/embed
+	cp Geo/GModel.cpp Geo/GModelIO_Mesh.cpp\
+             Geo/GEntity.cpp Geo/GVertex.cpp Geo/GEdge.cpp\
+             Geo/GEdgeLoop.cpp Geo/GFace.cpp Geo/GRegion.cpp\
+           Geo/MElement.cpp Geo/MFace.cpp Geo/MVertex.cpp\
+           Numeric/NumericEmbedded.cpp utils/embed
+
 variables: configure
 	@echo "********************************************************************"
 	@echo "Please configure Gmsh by running ./configure"
diff --git a/Numeric/Makefile b/Numeric/Makefile
index 5d40e84c35dc866d2a9584241e71f5d35436de81..f7398005017315de68d95f5012b7f2e07472894e 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.42 2008-01-19 22:06:03 geuzaine Exp $
+# $Id: Makefile,v 1.43 2008-01-22 16:47:10 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -26,6 +26,7 @@ INCLUDE = ${DASH}I../Common ${DASH}I../DataStr ${DASH}I../Numeric ${DASH}I../con
 CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE} ${SYSINCLUDE}
 
 SRC = Numeric.cpp\
+        NumericEmbedded.cpp\
       EigSolve.cpp\
       predicates.cpp\
       gsl_newt.cpp\
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index e32ee9249b434ac7dc76de761dbebc53c0f836ec..9b4341aa3bb035d7d42d588d55509f07ab9116a9 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1,4 +1,4 @@
-// $Id: Numeric.cpp,v 1.38 2008-01-19 22:06:04 geuzaine Exp $
+// $Id: Numeric.cpp,v 1.39 2008-01-22 16:47:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -19,9 +19,6 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
-// this file should contain only purely numerical routines (that do
-// not depend on any Gmsh structures)
-
 #include "Message.h"
 #include "Numeric.h"
 
@@ -80,215 +77,6 @@ int check_gsl()
 
 #endif
 
-double myatan2(double a, double b)
-{
-  if(a == 0.0 && b == 0)
-    return 0.0;
-  return atan2(a, b);
-}
-
-double myasin(double a)
-{
-  if(a <= -1.)
-    return -Pi / 2.;
-  else if(a >= 1.)
-    return Pi / 2.;
-  else
-    return asin(a);
-}
-
-double myacos(double a)
-{
-  if(a <= -1.)
-    return Pi;
-  else if(a >= 1.)
-    return 0.;
-  else
-    return acos(a);
-}
-
-
-void matvec(double mat[3][3], double vec[3], double res[3])
-{
-  res[0] = mat[0][0] * vec[0] + mat[0][1] * vec[1] + mat[0][2] * vec[2];
-  res[1] = mat[1][0] * vec[0] + mat[1][1] * vec[1] + mat[1][2] * vec[2];
-  res[2] = mat[2][0] * vec[0] + mat[2][1] * vec[1] + mat[2][2] * vec[2];
-}
-
-
-void normal3points(double x0, double y0, double z0,
-		   double x1, double y1, double z1,
-		   double x2, double y2, double z2,
-		   double n[3])
-{
-  double t1[3] = {x1 - x0, y1 - y0, z1 - z0};
-  double t2[3] = {x2 - x0, y2 - y0, z2 - z0};
-  prodve(t1, t2, n);
-  norme(n);
-}
-
-int sys2x2(double mat[2][2], double b[2], double res[2])
-{
-  double det, ud, norm;
-  int i;
-
-  norm = DSQR(mat[0][0]) + DSQR(mat[1][1]) + DSQR(mat[0][1]) + DSQR(mat[1][0]);
-  det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
-
-
-  // TOLERANCE ! WARNING WARNING
-  if(norm == 0.0 || fabs(det) / norm < 1.e-12) {
-    if(norm)
-      Msg(DEBUG, "Assuming 2x2 matrix is singular (det/norm == %.16g)",
-          fabs(det) / norm);
-    res[0] = res[1] = 0.0;
-    return 0;
-  }
-  ud = 1. / det;
-
-  res[0] = b[0] * mat[1][1] - mat[0][1] * b[1];
-  res[1] = mat[0][0] * b[1] - mat[1][0] * b[0];
-
-  for(i = 0; i < 2; i++)
-    res[i] *= ud;
-
-  return (1);
-}
-
-double det3x3(double mat[3][3])
-{
-  return (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
-	  mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
-	  mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]));
-}
-
-double trace3x3(double mat[3][3])
-{
-  return mat[0][0] + mat[1][1] + mat[2][2];
-}
-
-double trace2 (double mat[3][3])
-{
-  double a00 =  mat[0][0] * mat[0][0] + mat[1][0] * mat[0][1] + mat[2][0] * mat[0][2]; 
-  double a11 =  mat[1][0] * mat[0][1] + mat[1][1] * mat[1][1] + mat[1][2] * mat[2][1]; 
-  double a22 =  mat[2][0] * mat[0][2] + mat[2][1] * mat[1][2] + mat[2][2] * mat[2][2];
-  
-  return a00 + a11 + a22;
-}
-
-int sys3x3(double mat[3][3], double b[3], double res[3], double *det)
-{
-  double ud;
-  int i;
-
-  *det = det3x3(mat);
-
-  if(*det == 0.0) {
-    res[0] = res[1] = res[2] = 0.0;
-    return (0);
-  }
-
-  ud = 1. / (*det);
-
-  res[0] = b[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
-    mat[0][1] * (b[1] * mat[2][2] - mat[1][2] * b[2]) +
-    mat[0][2] * (b[1] * mat[2][1] - mat[1][1] * b[2]);
-
-  res[1] = mat[0][0] * (b[1] * mat[2][2] - mat[1][2] * b[2]) -
-    b[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
-    mat[0][2] * (mat[1][0] * b[2] - b[1] * mat[2][0]);
-
-  res[2] = mat[0][0] * (mat[1][1] * b[2] - b[1] * mat[2][1]) -
-    mat[0][1] * (mat[1][0] * b[2] - b[1] * mat[2][0]) +
-    b[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
-
-  for(i = 0; i < 3; i++)
-    res[i] *= ud;
-  return (1);
-
-}
-
-int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
-{
-  int out;
-  double norm;
-
-  out = sys3x3(mat, b, res, det);
-  norm =
-    DSQR(mat[0][0]) + DSQR(mat[0][1]) + DSQR(mat[0][2]) +
-    DSQR(mat[1][0]) + DSQR(mat[1][1]) + DSQR(mat[1][2]) +
-    DSQR(mat[2][0]) + DSQR(mat[2][1]) + DSQR(mat[2][2]);
-
-  // TOLERANCE ! WARNING WARNING
-  if(norm == 0.0 || fabs(*det) / norm < 1.e-12) {
-    if(norm)
-      Msg(DEBUG, "Assuming 3x3 matrix is singular (det/norm == %.16g)",
-          fabs(*det) / norm);
-    res[0] = res[1] = res[2] = 0.0;
-    return 0;
-  }
-
-  return out;
-}
-
-double det2x2(double mat[2][2])
-{
-  return mat[0][0] *mat[1][1] -mat[1][0] *mat[0][1];
-}
-
-double det2x3(double mat[2][3])
-{
-  double v1[3] = {mat[0][0],mat[0][1],mat[0][2]};
-  double v2[3] = {mat[1][0],mat[1][1],mat[1][2]};
-  double n[3];
-  
-  prodve (v1,v2,n);
-  return norm3(n);
-}
-
-double inv2x2(double mat[2][2], double inv[2][2])
-{
-  const double det = det2x2 ( mat );
-  if(det){
-    double ud = 1. / det;
-    inv[0][0] =  ( mat[1][1] ) * ud ;
-    inv[1][0] = -( mat[1][0] ) * ud ;
-    inv[0][1] = -( mat[0][1] ) * ud ;
-    inv[1][1] =  ( mat[0][0] ) * ud ;
-  }
-  else{
-    Msg(GERROR, "Singular matrix");
-    for(int i = 0; i < 2; i++)
-      for(int j = 0; j < 2; j++)
-	inv[i][j] = 0.;
-  }
-  return det;
-}
-
-double inv3x3(double mat[3][3], double inv[3][3])
-{
-  double det = det3x3(mat);
-  if(det){
-    double ud = 1. / det;
-    inv[0][0] =  ( mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1] ) * ud ;
-    inv[1][0] = -( mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0] ) * ud ;
-    inv[2][0] =  ( mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0] ) * ud ;
-    inv[0][1] = -( mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1] ) * ud ;
-    inv[1][1] =  ( mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0] ) * ud ;
-    inv[2][1] = -( mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0] ) * ud ;
-    inv[0][2] =  ( mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1] ) * ud ;
-    inv[1][2] = -( mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0] ) * ud ;
-    inv[2][2] =  ( mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0] ) * ud ;
-  }
-  else{
-    Msg(GERROR, "Singular matrix");
-    for(int i = 0; i < 3; i++)
-      for(int j = 0; j < 3; j++)
-	inv[i][j] = 0.;
-  }
-  return det;
-}
-
 void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
 {
   int i, j, k, n = 3;
@@ -361,236 +149,3 @@ void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
   free_dvector(W, 1, n);
 #endif
 }
-
-double angle_02pi(double A3)
-{
-  double DP = 2 * Pi;
-  while(A3 > DP || A3 < 0.) {
-    if(A3 > 0)
-      A3 -= DP;
-    else
-      A3 += DP;
-  }
-  return A3;
-}
-
-double angle_plan(double V[3], double P1[3], double P2[3], double n[3])
-{
-  double PA[3], PB[3], angplan;
-  double cosc, sinc, c[3];
-
-  PA[0] = P1[0] - V[0];
-  PA[1] = P1[1] - V[1];
-  PA[2] = P1[2] - V[2];
-
-  PB[0] = P2[0] - V[0];
-  PB[1] = P2[1] - V[1];
-  PB[2] = P2[2] - V[2];
-
-  norme(PA);
-  norme(PB);
-
-  prodve(PA, PB, c);
-
-  prosca(PA, PB, &cosc);
-  prosca(c, n, &sinc);
-  angplan = myatan2(sinc, cosc);
-
-  return angplan;
-}
-
-double triangle_area(double p0[3], double p1[3], double p2[3])
-{
-  double a[3], b[3], c[3];
-  
-  a[0] = p2[0] - p1[0];
-  a[1] = p2[1] - p1[1];
-  a[2] = p2[2] - p1[2];
-  
-  b[0] = p0[0] - p1[0];
-  b[1] = p0[1] - p1[1];
-  b[2] = p0[2] - p1[2];
-  
-  prodve(a, b, c);
-  return (0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]));
-}
-
-char float2char(float f)
-{
-  // float normalized in [-1, 1], char in [-127, 127]
-  f *= 127.;
-  if(f > 127.) return 127;
-  else if(f < -127.) return -127;
-  else return (char)f;
-}
-
-float char2float(char c)
-{
-  float f = c;
-  f /= 127.;
-  if(f > 1.) return 1.;
-  else if(f < -1) return -1.;
-  else return f;
-}
-
-double InterpolateIso(double *X, double *Y, double *Z,
-                      double *Val, double V, int I1, int I2,
-                      double *XI, double *YI, double *ZI)
-{
-  if(Val[I1] == Val[I2]) {
-    *XI = X[I1];
-    *YI = Y[I1];
-    *ZI = Z[I1];
-    return 0;
-  }
-  else {
-    double coef = (V - Val[I1]) / (Val[I2] - Val[I1]);
-    *XI = coef * (X[I2] - X[I1]) + X[I1];
-    *YI = coef * (Y[I2] - Y[I1]) + Y[I1];
-    *ZI = coef * (Z[I2] - Z[I1]) + Z[I1];
-    return coef;
-  }
-}
-
-void gradSimplex(double *x, double *y, double *z, double *v, double *grad)
-{
-  // p = p1 * (1-u-v-w) + p2 u + p3 v + p4 w
-
-  double mat[3][3];
-  double det, b[3];
-  mat[0][0] = x[1] - x[0];
-  mat[1][0] = x[2] - x[0];
-  mat[2][0] = x[3] - x[0];
-  mat[0][1] = y[1] - y[0];
-  mat[1][1] = y[2] - y[0];
-  mat[2][1] = y[3] - y[0];
-  mat[0][2] = z[1] - z[0];
-  mat[1][2] = z[2] - z[0];
-  mat[2][2] = z[3] - z[0];
-  b[0] = v[1] - v[0];
-  b[1] = v[2] - v[0];
-  b[2] = v[3] - v[0];
-  sys3x3(mat, b, grad, &det);
-}
-
-double ComputeVonMises(double *V)
-{
-  double tr = (V[0] + V[4] + V[8]) / 3.;
-  double v11 = V[0] - tr, v12 = V[1],      v13 = V[2];
-  double v21 = V[3],      v22 = V[4] - tr, v23 = V[5];
-  double v31 = V[6],      v32 = V[7],      v33 = V[8] - tr;
-  return sqrt(1.5 * (v11 * v11 + v12 * v12 + v13 * v13 +
-                     v21 * v21 + v22 * v22 + v23 * v23 +
-                     v31 * v31 + v32 * v32 + v33 * v33));
-}
-
-double ComputeScalarRep(int numComp, double *val)
-{
-  if(numComp == 1)
-    return val[0];
-  else if(numComp == 3)
-    return sqrt(val[0] * val[0] + val[1] * val[1] + val[2] * val[2]);
-  else if(numComp == 9)
-    return ComputeVonMises(val);
-  return 0.;
-}
-
-void eigenvalue(double mat[3][3], double v[3])
-{            
-  // characteristic polynomial of T : find v root of
-  // v^3 - I1 v^2 + I2 T - I3 = 0
-  // I1 : first invariant , trace(T)
-  // I2 : second invariant , 1/2 (I1^2 -trace(T^2))
-  // I3 : third invariant , det T
-
-  double c[4];
-  c[3] = 1.0;
-  c[2] = - trace3x3(mat);
-  c[1] = 0.5 * (c[2]*c[2] - trace2(mat));
-  c[0] = - det3x3(mat);
-  
-  //printf("%g %g %g\n", mat[0][0], mat[0][1], mat[0][2]);
-  //printf("%g %g %g\n", mat[1][0], mat[1][1], mat[1][2]);
-  //printf("%g %g %g\n", mat[2][0], mat[2][1], mat[2][2]);
-  //printf("%g x^3 + %g x^2 + %g x + %g = 0\n", c[3], c[2], c[1], c[0]);
-  
-  double imag[3];
-  FindCubicRoots(c, v, imag);
-  eigsort(v);
-}
-
-void FindCubicRoots(const double coef[4], double real[3], double imag[3])
-{
-  double a = coef[3];
-  double b = coef[2];
-  double c = coef[1];
-  double d = coef[0];
-
-  if(!a || !d){
-    Msg(GERROR, "Degenerate cubic: use a second degree solver!");
-    return;
-  }
-
-  b /= a;
-  c /= a;
-  d /= a;
-  
-  double q = (3.0*c - (b*b))/9.0;
-  double r = -(27.0*d) + b*(9.0*c - 2.0*(b*b));
-  r /= 54.0;
-
-  double discrim = q*q*q + r*r;
-  imag[0] = 0.0; // The first root is always real.
-  double term1 = (b/3.0);
-
-  if (discrim > 0) { // one root is real, two are complex
-    double s = r + sqrt(discrim);
-    s = ((s < 0) ? -pow(-s, (1.0/3.0)) : pow(s, (1.0/3.0)));
-    double t = r - sqrt(discrim);
-    t = ((t < 0) ? -pow(-t, (1.0/3.0)) : pow(t, (1.0/3.0)));
-    real[0] = -term1 + s + t;
-    term1 += (s + t)/2.0;
-    real[1] = real[2] = -term1;
-    term1 = sqrt(3.0)*(-t + s)/2;
-    imag[1] = term1;
-    imag[2] = -term1;
-    return;
-  }
-
-  // The remaining options are all real
-  imag[1] = imag[2] = 0.0;
-  
-  double r13;
-  if (discrim == 0){ // All roots real, at least two are equal.
-    r13 = ((r < 0) ? -pow(-r,(1.0/3.0)) : pow(r,(1.0/3.0)));
-    real[0] = -term1 + 2.0*r13;
-    real[1] = real[2] = -(r13 + term1);
-    return;
-  }
-
-  // Only option left is that all roots are real and unequal (to get
-  // here, q < 0)
-  q = -q;
-  double dum1 = q*q*q;
-  dum1 = acos(r/sqrt(dum1));
-  r13 = 2.0*sqrt(q);
-  real[0] = -term1 + r13*cos(dum1/3.0);
-  real[1] = -term1 + r13*cos((dum1 + 2.0*M_PI)/3.0);
-  real[2] = -term1 + r13*cos((dum1 + 4.0*M_PI)/3.0);
-}
-
-void  eigsort(double d[3])
-{
-  int k, j, i;
-  double p;
-  
-  for (i=0; i<3; i++) {
-    p=d[k=i];
-    for (j=i+1;j<3;j++)
-      if (d[j] >= p) p=d[k=j];
-    if (k != i) {
-      d[k]=d[i];
-      d[i]=p;
-    }
-  }
-}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 4051998d450ced4e5dd3e6ab49965c20ee79ecf0..5f926aceec060b0e36190aeb3383cbdc6475b6d7 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -20,127 +20,35 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
-#include <math.h>
-
-#define RADTODEG      57.295779513082321
-#define RacineDeDeux  1.4142135623730950
-#define RacineDeTrois 1.7320508075688773
-#define Pi            3.1415926535897932
-#define Deux_Pi       6.2831853071795865
-
-#if !defined(M_PI)
-#define M_PI Pi
-#endif
-
-#define MIN(a,b) (((a)<(b))?(a):(b))
-#define MAX(a,b) (((a)<(b))?(b):(a))
-#define SQR(a)   ((a)*(a))
-
-#define IMIN MIN
-#define LMIN MIN
-#define FMIN MIN
-#define DMIN MIN
-
-#define IMAX MAX
-#define LMAX MAX
-#define FMAX MAX
-#define DMAX MAX
-
-#define DSQR SQR
-#define FSQR SQR
-
-#define THRESHOLD(a,b,c) (((a)>(c))?(c):((a)<(b)?(b):(a)))
-
-#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
-#define sign(x)      (((x)>=0)?1:-1)
-#define Pred(x)      ((x)->prev)
-#define Succ(x)      ((x)->next)
-#define square(x)    ((x)*(x))
+#include "NumericEmbedded.h"
 
 int check_gsl();
-
-double myatan2(double a, double b);
-double myasin(double a);
-double myacos(double a);
-inline void prodve(double a[3], double b[3], double c[3])
-{
-  c[2] = a[0] * b[1] - a[1] * b[0];
-  c[1] = -a[0] * b[2] + a[2] * b[0];
-  c[0] = a[1] * b[2] - a[2] * b[1];
-}
-inline void prosca(double a[3], double b[3], double *c)
-{
-  *c = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
-}
-void matvec(double mat[3][3], double vec[3], double res[3]);
-inline double norm3(double a[3]){  return sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);}
-inline double norme(double a[3])
-{
-  const double mod = norm3(a);
-  if(mod != 0.0){
-    const double one_over_mod = 1./mod;
-    a[0] *= one_over_mod;
-    a[1] *= one_over_mod;
-    a[2] *= one_over_mod;
-  }
-  return mod;
-}
-
-void normal3points(double x0, double y0, double z0,
-		   double x1, double y1, double z1,
-		   double x2, double y2, double z2,
-		   double n[3]);
-int sys2x2(double mat[2][2], double b[2], double res[2]);
-int sys3x3(double mat[3][3], double b[3], double res[3], double *det);
-int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det);
-double det2x2(double mat[2][2]);
-double det2x3(double mat[2][3]);
-double det3x3(double mat[3][3]);
-double trace3x3(double mat[3][3]);
-double trace2 (double mat[3][3]);
-double inv3x3(double mat[3][3], double inv[3][3]);
-double inv2x2(double mat[2][2], double inv[2][2]);
 void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
-double angle_02pi(double A3);
-double angle_plan(double V[3], double P1[3], double P2[3], double n[3]);
-double triangle_area(double p0[3], double p1[3], double p2[3]);
-char float2char(float f);
-float char2float(char c);
-void eigenvalue(double mat[3][3], double re[3]);
-void FindCubicRoots(const double coeff[4], double re[3], double im[3]);
-void eigsort(double d[3]);
-
-double InterpolateIso(double *X, double *Y, double *Z, 
-		      double *Val, double V, int I1, int I2, 
-		      double *XI, double *YI ,double *ZI);
-void gradSimplex(double *x, double *y, double *z, double *v, double *grad);
-double ComputeVonMises(double *val);
-double ComputeScalarRep(int numComp, double *val);
-
-/* Numerical routines implemented using either Numerical Recipes or
-   the GSL */
 
+// Numerical routines implemented using either Numerical Recipes or
+// the GSL
 double brent(double ax, double bx, double cx,
 	     double (*f)(double), double tol, double *xmin);
 void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
 	    double *fc, double (*func)(double));
 void newt(double x[], int n, int *check,
 	  void (*vecfunc)(int, double [], double []));
-void minimize_2 ( double (*f) (double, double, void *data), 
-		  void (*df) (double, double, double &, double &, double &, void *data) ,
-		  void *data,int niter,
-		  double &U, double &V, double &res);
-void minimize_3 ( double (*f) (double, double, double, void *data), 
-		  void (*df) (double, double, double , double &, double &, double &, double &, void *data) ,
-		  void *data,int niter,
-		  double &U, double &V, double &W, double &res);
+void minimize_2 (double (*f) (double, double, void *data), 
+		 void (*df) (double, double, double &, double &, double &, void *data) ,
+		 void *data,int niter,
+		 double &U, double &V, double &res);
+void minimize_3 (double (*f) (double, double, double, void *data), 
+		 void (*df) (double, double, double , double &, double &, 
+			     double &, double &, void *data),
+		 void *data,int niter,
+		 double &U, double &V, double &W, double &res);
 void minimize_N (int N, 
 		 double (*f) (double*, void *data), 
 		 void (*df) (double*, double*, double &, void *data) ,
 		 void *data,int niter,
 		 double *, double &res);
 
-/* Robust geometrical predicates */
+// Robust geometrical predicates
 namespace gmsh{
 double exactinit();
 double incircle(double *pa, double *pb, double *pc, double *pd);
diff --git a/Numeric/NumericEmbedded.cpp b/Numeric/NumericEmbedded.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..06c28e6bcaffbefe375165fd648cd71988d8bece
--- /dev/null
+++ b/Numeric/NumericEmbedded.cpp
@@ -0,0 +1,470 @@
+// $Id: NumericEmbedded.cpp,v 1.1 2008-01-22 16:47:10 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+// this file should contain only purely numerical routines (that do
+// not depend on any Gmsh structures or external functions, expect
+// Msg)
+
+#include "NumericEmbedded.h"
+
+#if defined(HAVE_GMSH_EMBEDDED)
+#include "GmshEmbedded.h"
+#else
+#include "Message.h"
+#endif
+
+double myatan2(double a, double b)
+{
+  if(a == 0.0 && b == 0)
+    return 0.0;
+  return atan2(a, b);
+}
+
+double myasin(double a)
+{
+  if(a <= -1.)
+    return -Pi / 2.;
+  else if(a >= 1.)
+    return Pi / 2.;
+  else
+    return asin(a);
+}
+
+double myacos(double a)
+{
+  if(a <= -1.)
+    return Pi;
+  else if(a >= 1.)
+    return 0.;
+  else
+    return acos(a);
+}
+
+void matvec(double mat[3][3], double vec[3], double res[3])
+{
+  res[0] = mat[0][0] * vec[0] + mat[0][1] * vec[1] + mat[0][2] * vec[2];
+  res[1] = mat[1][0] * vec[0] + mat[1][1] * vec[1] + mat[1][2] * vec[2];
+  res[2] = mat[2][0] * vec[0] + mat[2][1] * vec[1] + mat[2][2] * vec[2];
+}
+
+void normal3points(double x0, double y0, double z0,
+		   double x1, double y1, double z1,
+		   double x2, double y2, double z2,
+		   double n[3])
+{
+  double t1[3] = {x1 - x0, y1 - y0, z1 - z0};
+  double t2[3] = {x2 - x0, y2 - y0, z2 - z0};
+  prodve(t1, t2, n);
+  norme(n);
+}
+
+int sys2x2(double mat[2][2], double b[2], double res[2])
+{
+  double det, ud, norm;
+  int i;
+
+  norm = DSQR(mat[0][0]) + DSQR(mat[1][1]) + DSQR(mat[0][1]) + DSQR(mat[1][0]);
+  det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
+
+  // TOLERANCE ! WARNING WARNING
+  if(norm == 0.0 || fabs(det) / norm < 1.e-12) {
+    if(norm)
+      Msg(DEBUG, "Assuming 2x2 matrix is singular (det/norm == %.16g)",
+          fabs(det) / norm);
+    res[0] = res[1] = 0.0;
+    return 0;
+  }
+  ud = 1. / det;
+
+  res[0] = b[0] * mat[1][1] - mat[0][1] * b[1];
+  res[1] = mat[0][0] * b[1] - mat[1][0] * b[0];
+
+  for(i = 0; i < 2; i++)
+    res[i] *= ud;
+
+  return (1);
+}
+
+double det3x3(double mat[3][3])
+{
+  return (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
+	  mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
+	  mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]));
+}
+
+double trace3x3(double mat[3][3])
+{
+  return mat[0][0] + mat[1][1] + mat[2][2];
+}
+
+double trace2 (double mat[3][3])
+{
+  double a00 =  mat[0][0] * mat[0][0] + mat[1][0] * mat[0][1] + mat[2][0] * mat[0][2]; 
+  double a11 =  mat[1][0] * mat[0][1] + mat[1][1] * mat[1][1] + mat[1][2] * mat[2][1]; 
+  double a22 =  mat[2][0] * mat[0][2] + mat[2][1] * mat[1][2] + mat[2][2] * mat[2][2];
+  
+  return a00 + a11 + a22;
+}
+
+int sys3x3(double mat[3][3], double b[3], double res[3], double *det)
+{
+  double ud;
+  int i;
+
+  *det = det3x3(mat);
+
+  if(*det == 0.0) {
+    res[0] = res[1] = res[2] = 0.0;
+    return (0);
+  }
+
+  ud = 1. / (*det);
+
+  res[0] = b[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
+    mat[0][1] * (b[1] * mat[2][2] - mat[1][2] * b[2]) +
+    mat[0][2] * (b[1] * mat[2][1] - mat[1][1] * b[2]);
+
+  res[1] = mat[0][0] * (b[1] * mat[2][2] - mat[1][2] * b[2]) -
+    b[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
+    mat[0][2] * (mat[1][0] * b[2] - b[1] * mat[2][0]);
+
+  res[2] = mat[0][0] * (mat[1][1] * b[2] - b[1] * mat[2][1]) -
+    mat[0][1] * (mat[1][0] * b[2] - b[1] * mat[2][0]) +
+    b[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
+
+  for(i = 0; i < 3; i++)
+    res[i] *= ud;
+  return (1);
+}
+
+int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
+{
+  int out;
+  double norm;
+
+  out = sys3x3(mat, b, res, det);
+  norm =
+    DSQR(mat[0][0]) + DSQR(mat[0][1]) + DSQR(mat[0][2]) +
+    DSQR(mat[1][0]) + DSQR(mat[1][1]) + DSQR(mat[1][2]) +
+    DSQR(mat[2][0]) + DSQR(mat[2][1]) + DSQR(mat[2][2]);
+
+  // TOLERANCE ! WARNING WARNING
+  if(norm == 0.0 || fabs(*det) / norm < 1.e-12) {
+    if(norm)
+      Msg(DEBUG, "Assuming 3x3 matrix is singular (det/norm == %.16g)",
+          fabs(*det) / norm);
+    res[0] = res[1] = res[2] = 0.0;
+    return 0;
+  }
+
+  return out;
+}
+
+double det2x2(double mat[2][2])
+{
+  return mat[0][0] *mat[1][1] -mat[1][0] *mat[0][1];
+}
+
+double det2x3(double mat[2][3])
+{
+  double v1[3] = {mat[0][0],mat[0][1],mat[0][2]};
+  double v2[3] = {mat[1][0],mat[1][1],mat[1][2]};
+  double n[3];
+  
+  prodve (v1,v2,n);
+  return norm3(n);
+}
+
+double inv2x2(double mat[2][2], double inv[2][2])
+{
+  const double det = det2x2 ( mat );
+  if(det){
+    double ud = 1. / det;
+    inv[0][0] =  ( mat[1][1] ) * ud ;
+    inv[1][0] = -( mat[1][0] ) * ud ;
+    inv[0][1] = -( mat[0][1] ) * ud ;
+    inv[1][1] =  ( mat[0][0] ) * ud ;
+  }
+  else{
+    Msg(GERROR, "Singular matrix");
+    for(int i = 0; i < 2; i++)
+      for(int j = 0; j < 2; j++)
+	inv[i][j] = 0.;
+  }
+  return det;
+}
+
+double inv3x3(double mat[3][3], double inv[3][3])
+{
+  double det = det3x3(mat);
+  if(det){
+    double ud = 1. / det;
+    inv[0][0] =  ( mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1] ) * ud ;
+    inv[1][0] = -( mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0] ) * ud ;
+    inv[2][0] =  ( mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0] ) * ud ;
+    inv[0][1] = -( mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1] ) * ud ;
+    inv[1][1] =  ( mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0] ) * ud ;
+    inv[2][1] = -( mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0] ) * ud ;
+    inv[0][2] =  ( mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1] ) * ud ;
+    inv[1][2] = -( mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0] ) * ud ;
+    inv[2][2] =  ( mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0] ) * ud ;
+  }
+  else{
+    Msg(GERROR, "Singular matrix");
+    for(int i = 0; i < 3; i++)
+      for(int j = 0; j < 3; j++)
+	inv[i][j] = 0.;
+  }
+  return det;
+}
+
+double angle_02pi(double A3)
+{
+  double DP = 2 * Pi;
+  while(A3 > DP || A3 < 0.) {
+    if(A3 > 0)
+      A3 -= DP;
+    else
+      A3 += DP;
+  }
+  return A3;
+}
+
+double angle_plan(double V[3], double P1[3], double P2[3], double n[3])
+{
+  double PA[3], PB[3], angplan;
+  double cosc, sinc, c[3];
+
+  PA[0] = P1[0] - V[0];
+  PA[1] = P1[1] - V[1];
+  PA[2] = P1[2] - V[2];
+
+  PB[0] = P2[0] - V[0];
+  PB[1] = P2[1] - V[1];
+  PB[2] = P2[2] - V[2];
+
+  norme(PA);
+  norme(PB);
+
+  prodve(PA, PB, c);
+
+  prosca(PA, PB, &cosc);
+  prosca(c, n, &sinc);
+  angplan = myatan2(sinc, cosc);
+
+  return angplan;
+}
+
+double triangle_area(double p0[3], double p1[3], double p2[3])
+{
+  double a[3], b[3], c[3];
+  
+  a[0] = p2[0] - p1[0];
+  a[1] = p2[1] - p1[1];
+  a[2] = p2[2] - p1[2];
+  
+  b[0] = p0[0] - p1[0];
+  b[1] = p0[1] - p1[1];
+  b[2] = p0[2] - p1[2];
+  
+  prodve(a, b, c);
+  return (0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]));
+}
+
+char float2char(float f)
+{
+  // float normalized in [-1, 1], char in [-127, 127]
+  f *= 127.;
+  if(f > 127.) return 127;
+  else if(f < -127.) return -127;
+  else return (char)f;
+}
+
+float char2float(char c)
+{
+  float f = c;
+  f /= 127.;
+  if(f > 1.) return 1.;
+  else if(f < -1) return -1.;
+  else return f;
+}
+
+double InterpolateIso(double *X, double *Y, double *Z,
+                      double *Val, double V, int I1, int I2,
+                      double *XI, double *YI, double *ZI)
+{
+  if(Val[I1] == Val[I2]) {
+    *XI = X[I1];
+    *YI = Y[I1];
+    *ZI = Z[I1];
+    return 0;
+  }
+  else {
+    double coef = (V - Val[I1]) / (Val[I2] - Val[I1]);
+    *XI = coef * (X[I2] - X[I1]) + X[I1];
+    *YI = coef * (Y[I2] - Y[I1]) + Y[I1];
+    *ZI = coef * (Z[I2] - Z[I1]) + Z[I1];
+    return coef;
+  }
+}
+
+void gradSimplex(double *x, double *y, double *z, double *v, double *grad)
+{
+  // p = p1 * (1-u-v-w) + p2 u + p3 v + p4 w
+
+  double mat[3][3];
+  double det, b[3];
+  mat[0][0] = x[1] - x[0];
+  mat[1][0] = x[2] - x[0];
+  mat[2][0] = x[3] - x[0];
+  mat[0][1] = y[1] - y[0];
+  mat[1][1] = y[2] - y[0];
+  mat[2][1] = y[3] - y[0];
+  mat[0][2] = z[1] - z[0];
+  mat[1][2] = z[2] - z[0];
+  mat[2][2] = z[3] - z[0];
+  b[0] = v[1] - v[0];
+  b[1] = v[2] - v[0];
+  b[2] = v[3] - v[0];
+  sys3x3(mat, b, grad, &det);
+}
+
+double ComputeVonMises(double *V)
+{
+  double tr = (V[0] + V[4] + V[8]) / 3.;
+  double v11 = V[0] - tr, v12 = V[1],      v13 = V[2];
+  double v21 = V[3],      v22 = V[4] - tr, v23 = V[5];
+  double v31 = V[6],      v32 = V[7],      v33 = V[8] - tr;
+  return sqrt(1.5 * (v11 * v11 + v12 * v12 + v13 * v13 +
+                     v21 * v21 + v22 * v22 + v23 * v23 +
+                     v31 * v31 + v32 * v32 + v33 * v33));
+}
+
+double ComputeScalarRep(int numComp, double *val)
+{
+  if(numComp == 1)
+    return val[0];
+  else if(numComp == 3)
+    return sqrt(val[0] * val[0] + val[1] * val[1] + val[2] * val[2]);
+  else if(numComp == 9)
+    return ComputeVonMises(val);
+  return 0.;
+}
+
+void eigenvalue(double mat[3][3], double v[3])
+{            
+  // characteristic polynomial of T : find v root of
+  // v^3 - I1 v^2 + I2 T - I3 = 0
+  // I1 : first invariant , trace(T)
+  // I2 : second invariant , 1/2 (I1^2 -trace(T^2))
+  // I3 : third invariant , det T
+
+  double c[4];
+  c[3] = 1.0;
+  c[2] = - trace3x3(mat);
+  c[1] = 0.5 * (c[2]*c[2] - trace2(mat));
+  c[0] = - det3x3(mat);
+  
+  //printf("%g %g %g\n", mat[0][0], mat[0][1], mat[0][2]);
+  //printf("%g %g %g\n", mat[1][0], mat[1][1], mat[1][2]);
+  //printf("%g %g %g\n", mat[2][0], mat[2][1], mat[2][2]);
+  //printf("%g x^3 + %g x^2 + %g x + %g = 0\n", c[3], c[2], c[1], c[0]);
+  
+  double imag[3];
+  FindCubicRoots(c, v, imag);
+  eigsort(v);
+}
+
+void FindCubicRoots(const double coef[4], double real[3], double imag[3])
+{
+  double a = coef[3];
+  double b = coef[2];
+  double c = coef[1];
+  double d = coef[0];
+
+  if(!a || !d){
+    Msg(GERROR, "Degenerate cubic: use a second degree solver!");
+    return;
+  }
+
+  b /= a;
+  c /= a;
+  d /= a;
+  
+  double q = (3.0*c - (b*b))/9.0;
+  double r = -(27.0*d) + b*(9.0*c - 2.0*(b*b));
+  r /= 54.0;
+
+  double discrim = q*q*q + r*r;
+  imag[0] = 0.0; // The first root is always real.
+  double term1 = (b/3.0);
+
+  if (discrim > 0) { // one root is real, two are complex
+    double s = r + sqrt(discrim);
+    s = ((s < 0) ? -pow(-s, (1.0/3.0)) : pow(s, (1.0/3.0)));
+    double t = r - sqrt(discrim);
+    t = ((t < 0) ? -pow(-t, (1.0/3.0)) : pow(t, (1.0/3.0)));
+    real[0] = -term1 + s + t;
+    term1 += (s + t)/2.0;
+    real[1] = real[2] = -term1;
+    term1 = sqrt(3.0)*(-t + s)/2;
+    imag[1] = term1;
+    imag[2] = -term1;
+    return;
+  }
+
+  // The remaining options are all real
+  imag[1] = imag[2] = 0.0;
+  
+  double r13;
+  if (discrim == 0){ // All roots real, at least two are equal.
+    r13 = ((r < 0) ? -pow(-r,(1.0/3.0)) : pow(r,(1.0/3.0)));
+    real[0] = -term1 + 2.0*r13;
+    real[1] = real[2] = -(r13 + term1);
+    return;
+  }
+
+  // Only option left is that all roots are real and unequal (to get
+  // here, q < 0)
+  q = -q;
+  double dum1 = q*q*q;
+  dum1 = acos(r/sqrt(dum1));
+  r13 = 2.0*sqrt(q);
+  real[0] = -term1 + r13*cos(dum1/3.0);
+  real[1] = -term1 + r13*cos((dum1 + 2.0*M_PI)/3.0);
+  real[2] = -term1 + r13*cos((dum1 + 4.0*M_PI)/3.0);
+}
+
+void  eigsort(double d[3])
+{
+  int k, j, i;
+  double p;
+  
+  for (i=0; i<3; i++) {
+    p=d[k=i];
+    for (j=i+1;j<3;j++)
+      if (d[j] >= p) p=d[k=j];
+    if (k != i) {
+      d[k]=d[i];
+      d[i]=p;
+    }
+  }
+}
diff --git a/Numeric/NumericEmbedded.h b/Numeric/NumericEmbedded.h
new file mode 100644
index 0000000000000000000000000000000000000000..091af663b078bfbd81d619c615cc8ee70fc034bf
--- /dev/null
+++ b/Numeric/NumericEmbedded.h
@@ -0,0 +1,118 @@
+#ifndef _NUMERIC_EMBEDDED_H_
+#define _NUMERIC_EMBEDDED_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include <math.h>
+
+#define RADTODEG      57.295779513082321
+#define RacineDeDeux  1.4142135623730950
+#define RacineDeTrois 1.7320508075688773
+#define Pi            3.1415926535897932
+#define Deux_Pi       6.2831853071795865
+
+#if !defined(M_PI)
+#define M_PI Pi
+#endif
+
+#define MIN(a,b) (((a)<(b))?(a):(b))
+#define MAX(a,b) (((a)<(b))?(b):(a))
+#define SQR(a)   ((a)*(a))
+
+#define IMIN MIN
+#define LMIN MIN
+#define FMIN MIN
+#define DMIN MIN
+
+#define IMAX MAX
+#define LMAX MAX
+#define FMAX MAX
+#define DMAX MAX
+
+#define DSQR SQR
+#define FSQR SQR
+
+#define THRESHOLD(a,b,c) (((a)>(c))?(c):((a)<(b)?(b):(a)))
+
+#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
+#define sign(x)      (((x)>=0)?1:-1)
+#define Pred(x)      ((x)->prev)
+#define Succ(x)      ((x)->next)
+#define square(x)    ((x)*(x))
+
+double myatan2(double a, double b);
+double myasin(double a);
+double myacos(double a);
+inline void prodve(double a[3], double b[3], double c[3])
+{
+  c[2] = a[0] * b[1] - a[1] * b[0];
+  c[1] = -a[0] * b[2] + a[2] * b[0];
+  c[0] = a[1] * b[2] - a[2] * b[1];
+}
+inline void prosca(double a[3], double b[3], double *c)
+{
+  *c = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+void matvec(double mat[3][3], double vec[3], double res[3]);
+inline double norm3(double a[3])
+{
+  return sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
+}
+inline double norme(double a[3])
+{
+  const double mod = norm3(a);
+  if(mod != 0.0){
+    const double one_over_mod = 1./mod;
+    a[0] *= one_over_mod;
+    a[1] *= one_over_mod;
+    a[2] *= one_over_mod;
+  }
+  return mod;
+}
+void normal3points(double x0, double y0, double z0,
+		   double x1, double y1, double z1,
+		   double x2, double y2, double z2,
+		   double n[3]);
+int sys2x2(double mat[2][2], double b[2], double res[2]);
+int sys3x3(double mat[3][3], double b[3], double res[3], double *det);
+int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det);
+double det2x2(double mat[2][2]);
+double det2x3(double mat[2][3]);
+double det3x3(double mat[3][3]);
+double trace3x3(double mat[3][3]);
+double trace2 (double mat[3][3]);
+double inv3x3(double mat[3][3], double inv[3][3]);
+double inv2x2(double mat[2][2], double inv[2][2]);
+double angle_02pi(double A3);
+double angle_plan(double V[3], double P1[3], double P2[3], double n[3]);
+double triangle_area(double p0[3], double p1[3], double p2[3]);
+char float2char(float f);
+float char2float(char c);
+void eigenvalue(double mat[3][3], double re[3]);
+void FindCubicRoots(const double coeff[4], double re[3], double im[3]);
+void eigsort(double d[3]);
+double InterpolateIso(double *X, double *Y, double *Z, 
+		      double *Val, double V, int I1, int I2, 
+		      double *XI, double *YI ,double *ZI);
+void gradSimplex(double *x, double *y, double *z, double *v, double *grad);
+double ComputeVonMises(double *val);
+double ComputeScalarRep(int numComp, double *val);
+
+#endif
diff --git a/utils/embed/GmshEmbedded.cpp b/utils/embed/GmshEmbedded.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d73949c0f4456f59b775dbe180516906d5a4b59
--- /dev/null
+++ b/utils/embed/GmshEmbedded.cpp
@@ -0,0 +1,84 @@
+#include "GModel.h"
+#include "GmshEmbedded.h"
+
+Context_T CTX;
+
+void GModel::createGEOInternals(){}
+void GModel::deleteGEOInternals(){}
+void GModel::deleteOCCInternals(){}
+
+void Msg(int level, const char *fmt, ...)
+{
+  va_list args;
+  int abort = 0;
+
+  va_start(args, fmt);
+
+  switch (level) {
+
+  case PROGRESS:
+  case STATUS1N:
+  case STATUS2N:
+    break;
+
+  case DIRECT:
+    if(CTX.verbosity >= 2) {
+      vfprintf(stdout, fmt, args);
+      fprintf(stdout, "\n");
+    }
+    break;
+
+  case FATAL:
+  case FATAL3: abort = 1;
+  case FATAL1:
+  case FATAL2:
+    fprintf(stderr, FATAL_STR);
+    vfprintf(stderr, fmt, args);
+    fprintf(stderr, "\n");
+    break;
+
+  case GERROR:
+  case GERROR1:
+  case GERROR2:
+  case GERROR3:
+    fprintf(stderr, ERROR_STR);
+    vfprintf(stderr, fmt, args);
+    fprintf(stderr, "\n");
+    break;
+
+  case WARNING:
+  case WARNING1:
+  case WARNING2:
+  case WARNING3:
+    if(CTX.verbosity >= 1) {
+      fprintf(stderr, WARNING_STR);
+      vfprintf(stderr, fmt, args);
+      fprintf(stderr, "\n");
+    }
+    break;
+
+  case DEBUG:
+  case DEBUG1:
+  case DEBUG2:
+  case DEBUG3:
+    if(CTX.verbosity >= 4) {
+      fprintf(stderr, DEBUG_STR);
+      vfprintf(stderr, fmt, args);
+      fprintf(stderr, "\n");
+    }
+    break;
+
+  default:
+    if(CTX.verbosity >= 2) {
+      fprintf(stderr, INFO_STR);
+      vfprintf(stderr, fmt, args);
+      fprintf(stderr, "\n");
+    }
+    break;
+  }
+
+  va_end(args);
+
+  if(abort)
+    exit(1);
+}
diff --git a/utils/embed/GmshEmbedded.h b/utils/embed/GmshEmbedded.h
new file mode 100644
index 0000000000000000000000000000000000000000..89586113ee0ad7c79f94fb0595713b7652c1164d
--- /dev/null
+++ b/utils/embed/GmshEmbedded.h
@@ -0,0 +1,95 @@
+#ifndef _GMSH_EMBEDDED_H_
+#define _GMSH_EMBEDDED_H_
+
+#include <stdarg.h>
+#include "NumericEmbedded.h"
+
+#define FATAL          1  // Fatal error (causes Gmsh to exit)
+#define FATAL1         2  // First part of a multiline FATAL message 
+#define FATAL2         3  // Middle part of a multiline FATAL message
+#define FATAL3         4  // Last part of a multiline FATAL message  
+
+#define GERROR         5  // Error (but Gmsh can live with it)
+#define GERROR1        6  // First part of a multiline ERROR message 
+#define GERROR2        7  // Middle part of a multiline ERROR message
+#define GERROR3        8  // Last part of a multiline ERROR message  
+
+#define WARNING        9  // Warning
+#define WARNING1      10  // First part of a multiline WARNING message 
+#define WARNING2      11  // Middle part of a multiline WARNING message
+#define WARNING3      12  // Last part of a multiline WARNING message  
+
+#define INFO          13  // Long informations
+#define INFO1         14  // First part of a multiline INFO message 
+#define INFO2         15  // Middle part of a multiline INFO message
+#define INFO3         16  // Last part of a multiline INFO message  
+
+#define DEBUG         17  // Long debug information
+#define DEBUG1        18  // First part of a multiline DEBUG message 
+#define DEBUG2        19  // Middle part of a multiline DEBUG message
+#define DEBUG3        20  // Last part of a multiline DEBUG message  
+
+#define STATUS1       21  // left status bar
+#define STATUS2       22  // right status bar
+
+#define STATUS1N      24  // Same as STATUS1, but not going into the log file
+#define STATUS2N      25  // Same as STATUS2, but not going into the log file
+
+#define ONSCREEN      27  // Persistent on-screen message
+
+#define DIRECT        30  // Direct message (no special formatting)
+#define SOLVER        31  // Solver message
+#define SOLVERR       32  // Solver errors and warnings
+
+#define PROGRESS      40  // Progress indicator
+
+#define WHITE_STR          "        : "
+#define FATAL_STR          "Fatal   : "
+#define ERROR_STR          "Error   : "
+#define WARNING_STR        "Warning : "
+#define INFO_STR           "Info    : "
+#define DEBUG_STR          "Debug   : "
+#define STATUS_STR         "Info    : "
+
+void   Signal(int signum);
+void   Msg(int level, const char *fmt, ...);
+
+class Context_T{
+public:
+  Context_T()
+  {
+    lc = 1.;
+    verbosity = 4;
+    hide_unselected = 0;
+    geom.tolerance = 1.e-6;
+    mesh.reverse_all_normals = 1;
+    pick_elements = 0;
+  }
+  double lc;
+  int verbosity;
+  int hide_unselected;
+  int pick_elements;
+  struct{
+    double tolerance;
+  } geom;
+  struct{
+    int reverse_all_normals;
+  } mesh;
+  unsigned int PACK_COLOR(int,int,int,int){ return 0; }
+  int UNPACK_RED(unsigned int){ return 0; }
+  int UNPACK_GREEN(unsigned int){ return 0; }
+  int UNPACK_BLUE(unsigned int){ return 0; }
+  int UNPACK_ALPHA(unsigned int){ return 0; }
+};
+
+class smooth_normals{
+ public:
+  int dummy;
+};
+
+class VertexArray{
+ public:
+  int dummy;
+};
+
+#endif
diff --git a/utils/embed/Makefile b/utils/embed/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..4d34695e59c2255b0c6809521ad726c7bea56412
--- /dev/null
+++ b/utils/embed/Makefile
@@ -0,0 +1,54 @@
+include ../../variables
+
+OBJEXT=.o
+LIBEXT=.a
+DASH=-
+OPTIM=-02
+FLAGS=
+SYSINCLUDE=
+
+
+LIB     = libGmsh${LIBEXT}
+INCLUDE = 
+CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE} ${SYSINCLUDE} ${DASH}DHAVE_GMSH_EMBEDDED
+
+SRC = GModel.cpp\
+        GModelIO_Mesh.cpp\
+      GEntity.cpp\
+        GVertex.cpp GEdge.cpp GEdgeLoop.cpp GFace.cpp GRegion.cpp\
+      MElement.cpp\
+        MFace.cpp MVertex.cpp\
+      NumericEmbedded.cpp\
+      GmshEmbedded.cpp
+
+OBJ = ${SRC:.cpp=${OBJEXT}}
+
+.SUFFIXES: ${OBJEXT} .cpp
+
+${LIB}: ${OBJ} 
+	${AR} ${ARFLAGS}${LIB} ${OBJ}
+	${RANLIB} ${LIB}
+
+.cpp${OBJEXT}:
+	${CXX} ${CFLAGS} ${DASH}c $<
+
+clean:
+	rm -f *${OBJEXT}
+
+distclean:
+	rm -f *${OBJEXT}
+	mkdir tmp
+	mv GmshEmbedded.* tmp
+	rm -f *.cpp *.h
+	mv tmp/* .
+	rmdir tmp
+
+depend:
+	(sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
+	${CXX} -MM ${CFLAGS} ${SRC} \
+	) >Makefile.new
+	cp Makefile Makefile.bak
+	cp Makefile.new Makefile
+	rm -f Makefile.new
+
+# DO NOT DELETE THIS LINE