diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index b53bf878083ee08a05cfbae44f95aaed024a4daf..414ffb0feaf344ee4aa85a840b59daceeb8ccbb0 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -700,7 +700,7 @@ void GFaceCompound::parametrize(iterationStep step) const
   Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step); 
   
   gmshAssembler<double> myAssembler(_lsys);
-  gmshFunction<double> ONE(1.0);
+  simpleFunction<double> ONE(1.0);
 
   if(_type == UNITCIRCLE){
     // maps the boundary onto a circle
@@ -849,8 +849,8 @@ void GFaceCompound::parametrize_conformal() const
     }    
   }    
 
-  gmshFunction<double> ONE(1.0);
-  gmshFunction<double> MONE(-1.0 );
+  simpleFunction<double> ONE(1.0);
+  simpleFunction<double> MONE(-1.0 );
   gmshLaplaceTerm laplace1(model(), &ONE, 1);
   gmshLaplaceTerm laplace2(model(), &ONE, 2);
   gmshCrossConfTerm cross12(model(), &ONE, 1 , 2);
@@ -900,7 +900,7 @@ void GFaceCompound::compute_distance() const
   double L  = norm(SVector3(bbox.max(), bbox.min())); 
   //printf("L=%g \n", L);
   double mu = L/28;
-  gmshFunction<double>  DIFF(mu*mu);
+  simpleFunction<double>  DIFF(mu*mu);
   gmshAssembler<double> myAssembler(_lsys);
   gmshDistanceTerm distance(model(), &DIFF, 1);
 
diff --git a/Numeric/gmshConvexCombination.h b/Numeric/gmshConvexCombination.h
index bb946e41d3668204b20be2b7361c3fc2881c357b..d3fc41a9483ddda9b6bc7d84e36973d55ae8fcee 100644
--- a/Numeric/gmshConvexCombination.h
+++ b/Numeric/gmshConvexCombination.h
@@ -7,7 +7,7 @@
 #define _GMSH_CONVEX_H_
 
 #include "gmshTermOfFormulation.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "Gmsh.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -15,7 +15,7 @@
 
 class gmshConvexCombinationTerm : public gmshNodalFemTerm<double> {
  protected:
-  const gmshFunction<double> *_diffusivity;
+  const simpleFunction<double> *_diffusivity;
   const int _iField ;
  public:
   virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
@@ -27,7 +27,7 @@ class gmshConvexCombinationTerm : public gmshNodalFemTerm<double> {
     *iFieldR = _iField;
   }
  public:
-  gmshConvexCombinationTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0) : 
+  gmshConvexCombinationTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0) : 
     gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField){}
   virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
 };
diff --git a/Numeric/gmshCrossConf.h b/Numeric/gmshCrossConf.h
index 792aac05ec6b566315e2b76d5016e67ef519f88a..dfdbe4d47cda2a430d6d59671baa3ae5e0ff09cf 100644
--- a/Numeric/gmshCrossConf.h
+++ b/Numeric/gmshCrossConf.h
@@ -9,7 +9,7 @@
 #define _GMSH_CROSSCONF_H_
 
 #include "gmshTermOfFormulation.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "Gmsh.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -17,7 +17,7 @@
 
 class gmshCrossConfTerm : public gmshNodalFemTerm<double> {
  protected:
-  const gmshFunction<double> *_diffusivity;
+  const simpleFunction<double> *_diffusivity;
   const int _iField;
   int _iFieldB ;
  public:
@@ -36,7 +36,7 @@ class gmshCrossConfTerm : public gmshNodalFemTerm<double> {
     *iFieldR = _iFieldB;
   }
  public:
-  gmshCrossConfTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) : 
+  gmshCrossConfTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) : 
     gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField)
     {
       _iFieldB = (iFieldB==-1) ? _iField : iFieldB;
@@ -47,7 +47,7 @@ class gmshCrossConfTerm : public gmshNodalFemTerm<double> {
 class gmshCrossConfTerm2DParametric : gmshCrossConfTerm {
   GFace *_gf;
  public:
-  gmshCrossConfTerm2DParametric(GFace *gf, gmshFunction<double> *diffusivity, int iField = 0) : 
+  gmshCrossConfTerm2DParametric(GFace *gf, simpleFunction<double> *diffusivity, int iField = 0) : 
     gmshCrossConfTerm(gf->model(), diffusivity, iField), _gf(gf) {}
   virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
 };
diff --git a/Numeric/gmshDistance.h b/Numeric/gmshDistance.h
index 57952ea7d4afc20bfc3d1a541e1fd1dfd6c2e1cf..af55a6fc67b5f4caee612eebf2f6521cdfd09fc9 100644
--- a/Numeric/gmshDistance.h
+++ b/Numeric/gmshDistance.h
@@ -9,7 +9,7 @@
 #define _GMSH_DISTANCE_H_
 
 #include "gmshTermOfFormulation.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "Gmsh.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -17,7 +17,7 @@
 
 class gmshDistanceTerm : public gmshNodalFemTerm<double> {
  protected:
-  const gmshFunction<double> *_diffusivity;
+  const simpleFunction<double> *_diffusivity;
   const int _iField;
   int _iFieldB ;
  public:
@@ -36,7 +36,7 @@ class gmshDistanceTerm : public gmshNodalFemTerm<double> {
     *iFieldR = _iFieldB;
   }
  public:
-  gmshDistanceTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) : 
+  gmshDistanceTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) : 
     gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField)
     {
       _iFieldB = (iFieldB==-1) ? _iField : iFieldB;
diff --git a/Numeric/gmshHelmholtz.h b/Numeric/gmshHelmholtz.h
index 800f5416a784d8e7dafe8ed264e35a68cd622daf..f2335f3bfe46da71f93fd4239389112e762a5f07 100644
--- a/Numeric/gmshHelmholtz.h
+++ b/Numeric/gmshHelmholtz.h
@@ -10,7 +10,7 @@
 
 #include <complex>
 #include "gmshTermOfFormulation.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "Gmsh.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -18,7 +18,7 @@
 
 class gmshHelmholtzTerm : public gmshNodalFemTerm<std::complex<double> > {
  private:
-  const gmshFunction<std::complex<double> > *_waveNumber;
+  const simpleFunction<std::complex<double> > *_waveNumber;
   const int _iField ;
  protected:
   virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
@@ -30,7 +30,7 @@ class gmshHelmholtzTerm : public gmshNodalFemTerm<std::complex<double> > {
     *iFieldR = _iField;
   }
  public:
-  gmshHelmholtzTerm(GModel *gm, gmshFunction<std::complex<double> > *waveNumber, 
+  gmshHelmholtzTerm(GModel *gm, simpleFunction<std::complex<double> > *waveNumber, 
                     int iField = 0) 
   : gmshNodalFemTerm<std::complex<double> >(gm), _waveNumber(waveNumber), _iField(iField){}
   virtual void elementMatrix(MElement *e, fullMatrix<std::complex<double> > &m) const;
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index 45c84f429fecdae55f7dcc60391808593e3ce5bb..cee00d054d042cbe3b1d066d436397d535268d84 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -9,7 +9,7 @@
 #define _GMSH_LAPLACE_H_
 
 #include "gmshTermOfFormulation.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "Gmsh.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -17,7 +17,7 @@
 
 class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
  protected:
-  const gmshFunction<double> *_diffusivity;
+  const simpleFunction<double> *_diffusivity;
   const int _iField;
   int _iFieldB ;
  public:
@@ -36,7 +36,7 @@ class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
     *iFieldR = _iFieldB;
   }
  public:
-  gmshLaplaceTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) : 
+  gmshLaplaceTerm(GModel *gm, simpleFunction<double> *diffusivity, int iField = 0, int iFieldB=-1) : 
     gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField)
     {
       _iFieldB = (iFieldB==-1) ? _iField : iFieldB;
diff --git a/Numeric/gmshProjection.h b/Numeric/gmshProjection.h
index 71287a8b6e982e51d23e5b237e92368dac1f2733..f20ef2f5ea983576032653d3ab21a58abd4f5e94 100644
--- a/Numeric/gmshProjection.h
+++ b/Numeric/gmshProjection.h
@@ -9,7 +9,7 @@
 #define _GMSH_PROJECTION_H_
 
 #include "gmshTermOfFormulation.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "Gmsh.h"
 #include "GModel.h"
 #include "MElement.h"
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index e6b7870488c73df4d83bc5e3ea31481c690ab6ac..721413c5e36832aea52243f66342adc8723e7165 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -12,7 +12,7 @@
 #include <map>
 #include <vector>
 #include "fullMatrix.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "gmshAssembler.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -107,7 +107,7 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
                     int dim,
                     int comp,
                     int field,
-                    const gmshFunction<scalar> &e,
+                    const simpleFunction<scalar> &e,
                     gmshAssembler<scalar> &lsys)
   {
     std::vector<MVertex *> v;
@@ -120,7 +120,7 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
                   int dim,
                   int comp,
                   int field,
-                  const gmshFunction<scalar> &fct,
+                  const simpleFunction<scalar> &fct,
                   gmshAssembler<scalar> &lsys)
   {
     std::map<int, std::vector<GEntity*> > groups[4];
diff --git a/Numeric/gmshFunction.h b/Numeric/simpleFunction.h
similarity index 59%
rename from Numeric/gmshFunction.h
rename to Numeric/simpleFunction.h
index f8d4fa5797621e0a8e5ff6152a8425329508c32f..b24057aef39a46604ce574eb3cc6ef0525d274f8 100644
--- a/Numeric/gmshFunction.h
+++ b/Numeric/simpleFunction.h
@@ -1,19 +1,17 @@
-// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
-
 // Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
-#ifndef _GMSH_FUNCTION_H_
-#define _GMSH_FUNCTION_H_
+#ifndef _SIMPLE_FUNCTION_H_
+#define _SIMPLE_FUNCTION_H_
 
 template <class scalar>
-class gmshFunction {
+class simpleFunction {
   scalar _val;
  public :
-  gmshFunction(scalar val=0) : _val(val) {}
-  virtual ~gmshFunction(){}
+  simpleFunction(scalar val=0) : _val(val) {}
+  virtual ~simpleFunction(){}
   virtual scalar operator () (double x, double y, double z) const { return _val; }
 };
 
diff --git a/Plugin/FiniteElement.cpp b/Plugin/FiniteElement.cpp
index 0b675655161109671be2fa6306292d2478ddbffc..0801b246d269acce97dd05dc10989ef224ef89af 100644
--- a/Plugin/FiniteElement.cpp
+++ b/Plugin/FiniteElement.cpp
@@ -68,7 +68,7 @@ StringXString *GMSH_FiniteElementPlugin::getOptionStr(int iopt)
 }
 
 template<class scalar>
-class gmshMathEvalFunction : public gmshFunction<scalar> {
+class gmshMathEvalFunction : public simpleFunction<scalar> {
  private:
   std::string _str;
   mathEvaluator *_f;
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index bbfed490deadb8806665e741d06b536e24ea8125..1ad666a2a92a5c9b659d2439e850083f9a1305f8 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -151,7 +151,7 @@ void elasticitySolver::solve()
        it != edgeDisplacements.end(); ++it){
     elasticityTerm El(pModel, 1, 1, _tag);
     El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag, 
-                        gmshFunction<double>(it->second), *pAssembler);
+                        simpleFunction<double>(it->second), *pAssembler);
     printf("-- Fixing edge %3d comp %1d to %8.5f\n", 
            it->first.first, it->first.second, it->second);
   }
@@ -160,7 +160,7 @@ void elasticitySolver::solve()
        it != faceDisplacements.end(); ++it){
     elasticityTerm El(pModel, 1, 1, _tag);
     El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag, 
-                        gmshFunction<double>(it->second), *pAssembler);
+                        simpleFunction<double>(it->second), *pAssembler);
     printf("-- Fixing face %3d comp %1d to %8.5f\n", 
            it->first.first, it->first.second, it->second);
   }
@@ -203,9 +203,9 @@ void elasticitySolver::solve()
     elasticityTerm El(pModel, 1, 1, _tag);
     int iEdge = it->first;
     SVector3 f = it->second;
-    El.neumannNodalBC(iEdge, 1, 0, _tag, gmshFunction<double>(f.x()), *pAssembler);
-    El.neumannNodalBC(iEdge, 1, 1, _tag, gmshFunction<double>(f.y()), *pAssembler);
-    El.neumannNodalBC(iEdge, 1, 2, _tag, gmshFunction<double>(f.z()), *pAssembler);
+    El.neumannNodalBC(iEdge, 1, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
+    El.neumannNodalBC(iEdge, 1, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
+    El.neumannNodalBC(iEdge, 1, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
     printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
   }
 
@@ -215,9 +215,9 @@ void elasticitySolver::solve()
     elasticityTerm El(pModel, 1, 1, _tag);
     int iFace = it->first;
     SVector3 f = it->second;
-    El.neumannNodalBC(iFace, 2, 0, _tag, gmshFunction<double>(f.x()), *pAssembler);
-    El.neumannNodalBC(iFace, 2, 1, _tag, gmshFunction<double>(f.y()), *pAssembler);
-    El.neumannNodalBC(iFace, 2, 2, _tag, gmshFunction<double>(f.z()), *pAssembler);
+    El.neumannNodalBC(iFace, 2, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
+    El.neumannNodalBC(iFace, 2, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
+    El.neumannNodalBC(iFace, 2, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
     printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
   }
   
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 51327f07048f1b0718764c079edc709fe994344c..4e042b3769b419c47c1b4f6c9f818b0cc35b3a02 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -10,7 +10,7 @@
 #include <map>
 #include <vector>
 #include "fullMatrix.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "dofManager.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -72,7 +72,7 @@ class femTerm {
     }
   }
   void dirichletNodalBC(int physical, int dim, int comp, int field,
-                        const gmshFunction<dataVec> &e,
+                        const simpleFunction<dataVec> &e,
                         dofManager<dataVec,dataMat> &dm)
   {
     std::vector<MVertex *> v;
@@ -82,7 +82,7 @@ class femTerm {
       dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
   }
   void neumannNodalBC(int physical, int dim, int comp,int field,
-                      const gmshFunction<dataVec> &fct,
+                      const simpleFunction<dataVec> &fct,
                       dofManager<dataVec,dataMat> &dm)
   {
     std::map<int, std::vector<GEntity*> > groups[4];
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 6475809459a1f134c6aec4f2bfc4ca1665f26328..1d094fede71f99b3d1683e971bd2c43f0d8699a9 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -8,7 +8,7 @@
 
 #include <assert.h>
 #include "femTerm.h"
-#include "gmshFunction.h"
+#include "simpleFunction.h"
 #include "Gmsh.h"
 #include "GModel.h"
 #include "MElement.h"
@@ -19,7 +19,7 @@
 template<class scalar>
 class helmoltzTerm : public femTerm<scalar, scalar> {
  protected:
-  const gmshFunction<scalar> *_k, *_a;
+  const simpleFunction<scalar> *_k, *_a;
   const int _iFieldR;
   int _iFieldC ;
  public:
@@ -36,8 +36,8 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
   }
   public:
   helmoltzTerm(GModel *gm, int iFieldR, int iFieldC, 
-               gmshFunction<scalar> *k, 
-               gmshFunction<scalar> *a) : 
+               simpleFunction<scalar> *k, 
+               simpleFunction<scalar> *a) : 
     femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR), _iFieldC(iFieldC){}
   virtual void elementMatrix(MElement *e, fullMatrix<scalar> &m) const
   {
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index 71b91f32a5ef375f9d42a78eea2ee0e534eefef9..f71288440198e520854b20a99b11a5a632367cdb 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -11,17 +11,17 @@
 // \nabla \cdot k \nabla U 
 template<class scalar> 
 class laplaceTerm : public helmholtzTerm<scalar> {
-  gmshFunction<scalar> *ONE;
+  simpleFunction<scalar> *ONE;
  public:
   laplaceTerm(GModel *gm, int iFieldR) : 
     helmholtzTerm<scalar>(gm, iFieldR, iFieldR, 0, 0)
   {
-    ONE = new  gmshFunction<scalar>(1.0);
+    ONE = new  simpleFunction<scalar>(1.0);
     _k = ONE;
   }
-  laplaceTerm(GModel *gm, int iFieldR, gmshFunction<scalar> *k) : 
+  laplaceTerm(GModel *gm, int iFieldR, simpleFunction<scalar> *k) : 
     helmholtzTerm<scalar>(gm, iFieldR, iFieldR, k, 0), ONE(0) {}
-  laplaceTerm(GModel *gm, int iFieldR, int iFieldC, gmshFunction<scalar> *k) : 
+  laplaceTerm(GModel *gm, int iFieldR, int iFieldC, simpleFunction<scalar> *k) : 
     helmholtzTerm<scalar>(gm, iFieldR, iFieldC, k, 0), ONE(0) {}
   virtual ~laplaceTerm()
   {