diff --git a/CMakeLists.txt b/CMakeLists.txt
index 92b4d421a9e3ff584b4187bd1aa1a408c041687e..6428f6ad3054af4550495762bcf85daeb5a56476 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -71,6 +71,7 @@ set(GMSH_API
     Mesh/meshGFaceDelaunayInsertion.h
   Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
     Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemFull.h 
+  Solver/elasticitySolver.h 
   Post/PView.h Post/PViewData.h Plugin/PluginManager.h
   Graphics/drawContext.h
   contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h 
@@ -487,11 +488,12 @@ if(ENABLE_MED OR ENABLE_CGNS)
 endif(ENABLE_MED OR ENABLE_CGNS)
 
 if(ENABLE_TAUCS)
-  find_library(TAUCS_LIB taucs)
+  find_library(TAUCS_LIB taucs PATHS ENV CASROOT PATH_SUFFIXES lib)
   if(TAUCS_LIB)
     find_path(TAUCS_INC "taucs.h" PATH_SUFFIXES src include)
     if(TAUCS_INC)
       set(HAVE_TAUCS TRUE)
+      add_definitions(-DTAUCS_CILK)
       list(APPEND CONFIG_OPTIONS "Taucs")
       list(APPEND EXTERNAL_LIBRARIES ${TAUCS_LIB})
       list(APPEND EXTERNAL_INCLUDES ${TAUCS_INC})
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 8532990ff5c07f6e04ccab57cefdef2db387e0dc..2fd1819ae0d00cc487e0d009a2090f1c2ebd60bf 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -331,13 +331,15 @@ void GFaceCompound::parametrize() const
     
     //Laplace parametrization
     //-----------------
-    parametrize(ITERU);
-    parametrize(ITERV);
-
+    if (_mapping == HARMONIC){
+      parametrize(ITERU);
+      parametrize(ITERV);
+    }
     //Conformal map parametrization
     //----------------- 
-    //parametrize_conformal();
-
+    else if (_mapping == CONFORMAL){
+      parametrize_conformal();
+    }
     //Distance function
     //-----------------
     //compute_distance();
@@ -476,9 +478,9 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
 GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
 			     std::list<GEdge*> &U0, std::list<GEdge*> &V0,
 			     std::list<GEdge*> &U1, std::list<GEdge*> &V1,
-			     linearSystem<double> *lsys)
+			     linearSystem<double> *lsys, typeOfMapping mpg)
   : GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0),
-    _lsys(lsys)
+    _lsys(lsys),_mapping(mpg)
 {
 
   if (!_lsys) {
@@ -925,7 +927,8 @@ void GFaceCompound::compute_distance() const
       SElement se((*it)->triangles[i]);
       distance.addToMatrix(myAssembler, &se);
     }
-    distance.addToRightHandSide(myAssembler, *it);
+    groupOfElements g(*it);
+    distance.addToRightHandSide(myAssembler, g);
   }
 
   Msg::Info("Distance Computation: Assembly done");
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 0a95efa9e3709bbefe553464f35466868fbb2f1b..6034b726137f30362a0d6cdf42bb29ed3c5fc7b7 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -46,6 +46,8 @@ class Octree;
 class GFaceCompound : public GFace {
  public:
   typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
+  typedef enum {HARMONIC=1,CONFORMAL=2} typeOfMapping;
+  typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
   void computeNormals(std::map<MVertex*, SVector3> &normals) const;
  protected:
   std::list<GFace*> _compound;
@@ -77,11 +79,11 @@ class GFaceCompound : public GFace {
   bool trivial() const ;
   linearSystem <double> *_lsys;
  public:
-  typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
   GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
                 std::list<GEdge*> &U0, std::list<GEdge*> &U1,
                 std::list<GEdge*> &V0, std::list<GEdge*> &V1,
-		linearSystem<double>* lsys =0);
+		linearSystem<double>* lsys =0,
+		typeOfMapping typ = HARMONIC);
   virtual ~GFaceCompound();
   Range<double> parBounds(int i) const 
   { return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
@@ -97,6 +99,7 @@ class GFaceCompound : public GFace {
   virtual int genusGeom ();
  private:
   typeOfIsomorphism _type;
+  typeOfMapping _mapping;
 };
 
 #endif
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 5937c3452d48e0487275c82c7d556728d3398095..b18781a80bee207d1d3fd26ff712b24e072865b6 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -172,8 +172,10 @@ int GModel::importGEOInternals()
       }
       GFace *gf = getFaceByTag(abs(p->Num));
       if (!gf){
-        GFaceCompound *gf = new GFaceCompound(this, p->Num, f_compound, 
-                                              b[0], b[1], b[2], b[3]);
+        GFaceCompound *gf = new GFaceCompound(this, abs(p->Num), f_compound, 
+                                              b[0], b[1], b[2], b[3],0,
+					      p->Num > 0 ? GFaceCompound::HARMONIC :
+					      GFaceCompound::CONFORMAL);
         add(gf);
       }
       else
diff --git a/Numeric/simpleFunction.h b/Numeric/simpleFunction.h
index b24057aef39a46604ce574eb3cc6ef0525d274f8..ecd2e520fe9f4678647b63401012147aad305e26 100644
--- a/Numeric/simpleFunction.h
+++ b/Numeric/simpleFunction.h
@@ -13,6 +13,9 @@ class simpleFunction {
   simpleFunction(scalar val=0) : _val(val) {}
   virtual ~simpleFunction(){}
   virtual scalar operator () (double x, double y, double z) const { return _val; }
+  virtual void gradient (double x, double y, double z,
+			 scalar & dfdx, scalar & dfdy, scalar & dfdz) const 
+  { dfdx = dfdy = dfdz = 0.0; }
 };
 
 #endif
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index a115f56ef491d2d411204de5fd3cb4f963ad0867..89621c741947718f1b7d5f75af4eede523b11d76 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -6,8 +6,10 @@
 set(SRC
   linearSystemCSR.cpp
   linearSystemTAUCS.cpp
+  groupOfElements.cpp
   elasticityTerm.cpp
   elasticitySolver.cpp
+  SElement.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 
diff --git a/Solver/SElement.h b/Solver/SElement.h
index b9a658f4bc7ff6c6e2a54b3445b4efdc08f9db46..518800fefac80ef95aa0130cefc28d50b477ec8c 100644
--- a/Solver/SElement.h
+++ b/Solver/SElement.h
@@ -7,6 +7,7 @@
 #define _SELEMENT_H_
 
 #include "MElement.h"
+#include "simpleFunction.h"
 
 // A solver element.
 
@@ -17,10 +18,29 @@ class SElement
   MElement *_e;
   // store discrete function space and other data here
   // ...
+  static simpleFunction<double> *_enrichement_s,*_enrichement_t;
+  // gradient of functions (possibly enriched)
+  void nodalFunctions ( double u, double v, double w, double s[],
+			simpleFunction<double> *_enrichement);   
+  void gradNodalFunctions ( double u, double v, double w, double invjac[3][3],
+			    double grad[][3], simpleFunction<double> *_enrichment);   
  public:
-  SElement(MElement *e) : _e(e) {}
+ SElement(MElement *e) : _e(e) {}  
   ~SElement(){}
-  MElement *getMeshElement() const { return _e; }
+  MElement *getMeshElement() const { return _e; }  
+  static void setShapeEnrichement(simpleFunction<double>*f){_enrichement_s=f;}
+  static void setTestEnrichement(simpleFunction<double>*f){_enrichement_t=f;}
+  static const simpleFunction<double>* getShapeEnrichement(){return _enrichement_s;}
+  static const simpleFunction<double>* getTestEnrichement(){return _enrichement_t;}
+  int getNumNodalShapeFunctions() const;
+  inline MVertex *getVertex(int i) const {return _e->getParent() ? _e->getParent()->getVertex(i) : _e->getVertex(i) ; }
+  int getNumNodalTestFunctions() const;
+  void nodalShapeFunctions ( double u, double v, double w, double s[]); 
+  void gradNodalShapeFunctions ( double u, double v, double w, double invjac[3][3],
+				 double grad[][3]); 
+  void nodalTestFunctions ( double u, double v, double w, double s[]); 
+  void gradNodalTestFunctions ( double u, double v, double w, double invjac[3][3],
+			        double grad[][3]);   
 };
 
 #endif
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 5410382b573630c832b57555f15bf0a72b7f148c..8ff1059539c02a96ea5f91463c73b98fe28d560f 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -4,12 +4,48 @@
 // bugs and problems to <gmsh@geuz.org>.
 
 #include <string.h>
+#include "GmshConfig.h"
 #include "elasticityTerm.h"
 #include "MTriangle.h"
 #include "MTetrahedron.h"
 #include "elasticitySolver.h"
 #include "linearSystemTAUCS.h"
+#include "linearSystemGMM.h"
 #include "Numeric.h"
+#include "PView.h"
+#include "PViewData.h"
+
+PView* elasticitySolver::buildDisplacementView  (const std::string &postFileName){
+  std::map<int, std::vector<double> > data;
+  std::set<MVertex*> v;
+  for (int i=0;i<elasticFields.size();++i){
+    groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+    for ( ; it != elasticFields[i].g->end() ; ++it ){
+      SElement se(*it);
+      for (int j=0;j<se.getNumNodalShapeFunctions();++j){
+	v.insert(se.getVertex(j));
+      }
+    }
+  }
+
+  std::set<MVertex*>::iterator it = v.begin();
+  for ( ; it!=v.end();++it){
+    std::vector<double> d;
+    d.push_back(0);d.push_back(0);d.push_back(0);
+    for (int i=0;i<elasticFields.size();++i){
+      const double E = (elasticFields[i]._enrichment) ? (*elasticFields[i]._enrichment)((*it)->x(),(*it)->y(),(*it)->z()) : 1.0;	
+      //      printf("%g %d \n",pAssembler->getDofValue(*it,0,elasticFields[i]._tag),elasticFields[i]._tag);
+      d[0] += pAssembler->getDofValue(*it,0,elasticFields[i]._tag)  * E;
+      d[1] += pAssembler->getDofValue(*it,1,elasticFields[i]._tag)  * E;
+      d[2] += pAssembler->getDofValue(*it,2,elasticFields[i]._tag)  * E;
+    }
+    data[(*it)->getNum()] = d;
+  }
+
+  PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
+  return pv;  
+}
+
 
 static double vonMises(dofManager<double,double> *a, MElement *e, 
                        double u, double v, double w, 
@@ -65,11 +101,21 @@ void elasticitySolver::readInputFile(const std::string &fn)
   while(!feof(f)){
     if(fscanf(f, "%s", what) != 1) return;
     // printf("%s\n", what);
-    if (!strcmp(what,"ElasticMaterial")){
-      double E, nu;
+    if (!strcmp(what,"ElasticDomain")){
+      elasticField field;
       int physical;
-      if(fscanf(f, "%d %lf %lf", &physical, &E, &nu) != 3) return;
-      elasticConstants[physical] = std::make_pair(E, nu);    
+      if(fscanf(f, "%d %lf %lf", &physical, &field._E, &field._nu) != 3) return;
+      field._tag = _tag;
+      //      printf("E %g nu %g\n",field._E,field._nu);
+      field.g = new groupOfElements (_dim, physical);
+      elasticFields.push_back(field);
+    }
+    else if (!strcmp(what,"Void")){
+      //      elasticField field;
+      //      create enrichment ...
+      //      create the group ...
+      //      assign a tag
+      //      elasticFields.push_back(field);
     }
     else if (!strcmp(what, "NodalDisplacement")){
       double val;
@@ -128,7 +174,8 @@ void elasticitySolver::solve()
 #ifdef HAVE_TAUCS
   linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
 #else
-  linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
+  linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+  lsys->setNoisy(2);
 #endif
   pAssembler = new dofManager<double, double>(lsys);
   
@@ -167,16 +214,15 @@ void elasticitySolver::solve()
   
   // we number the nodes : when a node is numbered, it cannot be numbered
   // again with another number. so we loop over elements
-  for (std::map<int, std::pair<double, double> >::iterator it = elasticConstants.begin();
-       it != elasticConstants.end() ; it++){
-    int physical = it->first;
-    std::vector<MVertex *> v;     
-    pModel->getMeshVerticesForPhysicalGroup(_dim, physical, v);
-    printf("Physical %d, dim: %d, nb vert: %d\n", physical, _dim, (int)v.size());
-    for (unsigned int i = 0; i < v.size(); i++){  
-      pAssembler->numberVertex(v[i], 0, _tag);  
-      pAssembler->numberVertex(v[i], 1, _tag);  
-      pAssembler->numberVertex(v[i], 2, _tag);  
+  for (int i=0;i<elasticFields.size();++i){
+    groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+    for ( ; it != elasticFields[i].g->end() ; ++it ){
+      SElement se(*it);
+      for (int j=0;j<se.getNumNodalShapeFunctions();++j){
+	pAssembler->numberVertex(se.getVertex(j), 0, elasticFields[i]._tag);  
+	pAssembler->numberVertex(se.getVertex(j), 1, elasticFields[i]._tag);  
+	pAssembler->numberVertex(se.getVertex(j), 2, elasticFields[i]._tag);  	
+      }
     }
   }
   
@@ -231,23 +277,28 @@ void elasticitySolver::solve()
     printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
     std::vector<GEntity*> ent = groups[_dim][iVolume];
     for (unsigned int i = 0; i < ent.size(); i++){
-      El.addToRightHandSide(*pAssembler, ent[i]);
+      //   to do
+      //      El.addToRightHandSide(*pAssembler, ent[i]);
     }
   }
   
   // assembling the stifness matrix
-  for (std::map<int, std::pair<double, double> >::iterator it = elasticConstants.begin();
-       it != elasticConstants.end() ; it++){
-    int physical = it->first;
-    double E = it->second.first;
-    double nu = it->second.second;
-    elasticityTerm El(0, E, nu, _tag);
-    std::vector<GEntity*> ent = groups[_dim][physical];
-    for (unsigned int i = 0; i < ent.size(); i++){
-      El.addToMatrix(*pAssembler, ent[i]);
+  for (int i=0;i<elasticFields.size();i++){
+    SElement::setShapeEnrichement(elasticFields[i]._enrichment);
+    for (int j=0;j<elasticFields.size();j++){
+      elasticityTerm El(0, elasticFields[i]._E, elasticFields[i]._nu,
+			elasticFields[i]._tag,elasticFields[j]._tag);
+      SElement::setTestEnrichement(elasticFields[j]._enrichment);
+      //      printf("%d %d\n",elasticFields[i]._tag,elasticFields[j]._tag);
+      El.addToMatrix(*pAssembler,*elasticFields[i].g, *elasticFields[j].g);      
     }
   }
 
+  //  for (int i=0;i<pAssembler->sizeOfR();i++){
+    //    printf("%g ",lsys->getFromRightHandSide(i));
+  //  }
+  //  printf("\n");
+
   // solving
   lsys->systemSolve();
 }
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index b59abbb565338e9e97e71abfef1ed6d386086c99..eba195515e69c6e8e6ba3af4c58451ffc9d93912 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -10,8 +10,19 @@
 #include <string>
 #include "SVector3.h"
 #include "dofManager.h"
+#include "simpleFunction.h"
 
 class GModel;
+class PView;
+class groupOfElements;
+
+struct elasticField {
+  int _tag; // tag for the dofManager
+  groupOfElements *g; // support for this field
+  simpleFunction<double> *_enrichment; // XFEM enrichment
+  double _E,_nu; // specific elastic datas (should be somewhere else)
+  elasticField () : g(0),_enrichment(0){}
+};
 
 // an elastic solver ...
 class elasticitySolver{
@@ -20,7 +31,7 @@ class elasticitySolver{
   int _dim, _tag;
   dofManager<double, double> *pAssembler;
   // young modulus and poisson coefficient per physical
-  std::map<int, std::pair<double, double> > elasticConstants;
+  std::vector<elasticField > elasticFields;
   // imposed nodal forces
   std::map<int, SVector3> nodalForces;
   // imposed line forces
@@ -45,14 +56,14 @@ class elasticitySolver{
   {
     nodalDisplacements[std::make_pair(iNode, dir)] = val;
   }
-  void addElasticConstants(double e, double nu, int physical)
-  {
-    elasticConstants[physical] = std::make_pair(e, nu);
-  }
+  //  void addElasticConstants(double e, double nu, int physical)
+  //  {
+  //    elasticConstants[physical] = std::make_pair(e, nu);
+  //  }
   void setMesh(const std::string &meshFileName);  
   virtual void solve();  
   void readInputFile(const std::string &meshFileName);
-  // PView *buildDisplacementView(const std::string &postFileName);
+  PView *buildDisplacementView(const std::string &postFileName);
   // PView *buildVonMisesView(const std::string &postFileName);
   // std::pair<PView *, PView*> buildErrorEstimateView
   //   (const std::string &errorFileName, double, int);
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index d67462bafd2d1ee25632edad7bc6cc43c4982f23..5b42a5227a7667634c56dbac381306dd64e2c955 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -6,16 +6,20 @@
 #include "elasticityTerm.h"
 #include "Numeric.h"
 
+// The SElement (Solver element) that has been sent to the function
+// contains 2 enrichments, that can enrich both shape and test functions   
+
 void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
 {
   MElement *e = se->getMeshElement();
-  int nbNodes = e->getNumVertices();
+  int nbNodes = se->getNumNodalShapeFunctions();
   int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
   int npts;
   IntPt *GP;
   double jac[3][3];
   double invjac[3][3];
-  double Grads[256][3], grads[100][3];
+  double Grads[100][3];
+  double GradsT[100][3];
   e->getIntegrationPoints(integrationOrder, &npts, &GP);
   
   m.set_all(0.);
@@ -46,34 +50,62 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
     const double w = GP[i].pt[2];
     const double weight = GP[i].weight;
     const double detJ = e->getJacobian(u, v, w, jac);
-    e->getGradShapeFunctions(u, v, w, grads);
     inv3x3(jac, invjac);
+    se->gradNodalShapeFunctions(u, v, w, invjac,Grads);
+
     B.set_all(0.);
     BT.set_all(0.);
-    for (int j = 0; j < nbNodes; j++){
-      Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] + 
-        invjac[0][2] * grads[j][2];
-      Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] + 
-        invjac[1][2] * grads[j][2];
-      Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] + 
-        invjac[2][2] * grads[j][2];
-      
-      BT(j, 0) = B(0, j) = Grads[j][0];
-      BT(j, 3) = B(3, j) = Grads[j][1];
-      BT(j, 4) = B(4, j) = Grads[j][2];
-      
-      BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
-      BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
-      BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
-      
-      BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
-      BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
-      BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
+    
+    if (se->getShapeEnrichement() == se->getTestEnrichement()){
+      for (int j = 0; j < nbNodes; j++){
+
+	//	printf(" GR(j) = %12.5E,%12.5E,%12.5E\n", Grads[j][0],Grads[j][1],Grads[j][2]);
+
+	BT(j, 0) = B(0, j) = Grads[j][0];
+	BT(j, 3) = B(3, j) = Grads[j][1];
+	BT(j, 4) = B(4, j) = Grads[j][2];
+	
+	BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
+	BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
+	BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
+	
+	BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
+	BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
+	BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
+      }
+    }
+    else{
+      printf("coucou AAAAAAAAAAAAARGH \n");
+      se->gradNodalTestFunctions (u, v, w, invjac,GradsT);
+      for (int j = 0; j < nbNodes; j++){
+	BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0];
+	BT(j, 3) = Grads[j][1]; B(3, j) = GradsT[j][1];
+	BT(j, 4) = Grads[j][2]; B(4, j) = GradsT[j][2];
+	
+	BT(j + nbNodes, 1) = Grads[j][1]; B(1, j + nbNodes) = GradsT[j][1];
+	BT(j + nbNodes, 3) = Grads[j][0]; B(3, j + nbNodes) = GradsT[j][0];
+	BT(j + nbNodes, 5) = Grads[j][2]; B(5, j + nbNodes) = GradsT[j][2];
+	
+	BT(j + 2 * nbNodes, 2) = Grads[j][2]; B(2, j + 2 * nbNodes) = GradsT[j][2];
+	BT(j + 2 * nbNodes, 4) = Grads[j][0]; B(4, j + 2 * nbNodes) = GradsT[j][0];
+	BT(j + 2 * nbNodes, 5) = Grads[j][1]; B(5, j + 2 * nbNodes) = GradsT[j][1];
+      }
     }
+
     BTH.set_all(0.);
     BTH.gemm(BT, H); 
     m.gemm(BTH, B, weight * detJ, 1.);
-  } 
+  }
+  return;
+  for (int i=0;i<3 * nbNodes;i++){
+    for (int j=0;j<3 * nbNodes;j++){
+      printf("%g ",m(i,j));
+    }
+    printf("\n");
+  }
+    printf("\n");
+    printf("\n");
+    
 }
 
 void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
@@ -95,7 +127,7 @@ void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
     const double w = GP[i].pt[2];
     const double weight = GP[i].weight;
     const double detJ = e->getJacobian(u, v, w, jac);
-    e->getShapeFunctions(u, v, w, ff);
+    se->nodalTestFunctions(u, v, w, ff);
     for (int j = 0; j < nbNodes; j++){
       m(j) += _volumeForce.x() *ff[j] * weight * detJ *.5;
       m(j + nbNodes) += _volumeForce.y() *ff[j] * weight * detJ *.5;
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index a03d459cdb9fc59cd1d263c18958df12d85e367a..ae826d930c3a9386ad6c6035686d9d5d12d94f54 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -15,9 +15,11 @@
 class elasticityTerm : public femTerm<double, double> {
  protected:
   double _E, _nu;
-  int _iField;
+  int _iFieldR, _iFieldC;
   SVector3 _volumeForce;
  public:
+  void setFieldC(int i){_iFieldC = i;} 
+  void setFieldR(int i){_iFieldR = i;} 
   // element matrix size : 3 dofs per vertex
   virtual int sizeOfR(SElement *se) const 
   {
@@ -35,11 +37,21 @@ class elasticityTerm : public femTerm<double, double> {
     int iCompR = iRow / e->getNumVertices();
     int ithLocalVertex = iRow % e->getNumVertices();
     return Dof(e->getVertex(ithLocalVertex)->getNum(),
-               Dof::createTypeWithTwoInts(iCompR, _iField));
+               Dof::createTypeWithTwoInts(iCompR, _iFieldR));
+  }
+  Dof getLocalDofC(SElement *se, int iCol) const 
+  {
+    MElement *e = se->getMeshElement();
+    int iCompC = iCol / e->getNumVertices();
+    int ithLocalVertex = iCol % e->getNumVertices();
+    return Dof(e->getVertex(ithLocalVertex)->getNum(),
+               Dof::createTypeWithTwoInts(iCompC, _iFieldC));
   }
  public:
-  elasticityTerm(GModel *gm, double E, double nu, int iField)
-    : femTerm<double, double>(gm), _E(E), _nu(nu), _iField(iField) {}
+ elasticityTerm(GModel *gm, double E, double nu, int fieldr, int fieldc)
+   : femTerm<double, double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldc) {}
+ elasticityTerm(GModel *gm, double E, double nu, int fieldr)
+   : femTerm<double, double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldr) {}
   void setVector(const SVector3 &f) { _volumeForce = f; }
   void elementMatrix(SElement *se, fullMatrix<double> &m) const;
   void elementVector(SElement *se, fullVector<double> &m) const;
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 29372dff928a9fca050181b7784b5d2bdb8a7abb..65c6938f5fa719d010b1c409b994ea7403ab7d92 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -14,6 +14,7 @@
 #include "dofManager.h"
 #include "GModel.h"
 #include "SElement.h"
+#include "groupOfElements.h"
 
 // a nodal finite element term : variables are always defined at nodes
 // of the mesh
@@ -42,15 +43,23 @@ class femTerm {
   {
     m.scale(0.0);
   }
-  // add the contribution from all the elements in the entity ge to
-  // the dof manager
-  void addToMatrix(dofManager<dataVec, dataMat> &dm, GEntity *ge) const
+
+  // add the contribution from all the elements in the intersection
+  // of two element groups L and C
+  void addToMatrix(dofManager<dataVec, dataMat> &dm, 
+		   groupOfElements &L, 
+		   groupOfElements &C) const
   {
-    for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
-      SElement se(ge->getMeshElement(i));
-      addToMatrix(dm, &se);
+    groupOfElements::elementContainer::const_iterator it = L.begin();
+    for ( ; it != L.end() ; ++it){
+      MElement *eL = *it;
+      if ( &C == &L || C.find(eL) ){
+	SElement se(eL);
+	addToMatrix(dm, &se);
+      }
     }
   }
+
   // add the contribution from a single element to the dof manager
   void addToMatrix(dofManager<dataVec, dataMat> &dm, SElement *se) const
   {
@@ -120,10 +129,13 @@ class femTerm {
       }
     }
   }
-  void addToRightHandSide(dofManager<dataVec, dataMat> &dm, GEntity *ge) const 
+
+  void addToRightHandSide(dofManager<dataVec, dataMat> &dm, groupOfElements &C) const 
   {
-    for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
-      SElement se(ge->getMeshElement(i));
+    groupOfElements::elementContainer::const_iterator it = C.begin();
+    for ( ; it != C.end() ; ++it){
+      MElement *eL = *it;
+      SElement se(eL);
       int nbR = sizeOfR(&se);
       fullVector<dataVec> V(nbR);
       elementVector(&se, V);
diff --git a/Solver/groupOfElements.cpp b/Solver/groupOfElements.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8fcaea61f529e070857efb9302aed6755f5ffd93
--- /dev/null
+++ b/Solver/groupOfElements.cpp
@@ -0,0 +1,29 @@
+#include "groupOfElements.h"
+#include "GModel.h"
+#include "GEntity.h"
+
+groupOfElements::groupOfElements(GFace*gf)
+{
+  elementFilterTrivial filter;
+  addElementary(gf,filter);  
+} 
+
+void groupOfElements::addElementary(GEntity *ge, const elementFilter &filter){
+  for (int j=0;j<ge->getNumMeshElements();j++){
+    MElement *e = ge->getMeshElement(j);
+    if (filter(e)){
+      insert(e);
+    }
+  }
+}
+
+void groupOfElements::addPhysical(int dim, int physical,
+				  const elementFilter &filter){
+  std::map<int, std::vector<GEntity*> > groups[4];
+  GModel::current()->getPhysicalGroups(groups);
+  std::vector<GEntity*> ent = groups[dim][physical];
+  for (int i=0;i<ent.size();i++){
+    addElementary(ent[i],filter);
+  }
+}
+
diff --git a/Solver/groupOfElements.h b/Solver/groupOfElements.h
new file mode 100644
index 0000000000000000000000000000000000000000..eae4e9eb58a11d3db9fefafc5e7916ebe84a54a8
--- /dev/null
+++ b/Solver/groupOfElements.h
@@ -0,0 +1,79 @@
+#ifndef _GROUPOFELEMENTS_H_
+#define _GROUPOFELEMENTS_H_
+
+#include <set>
+#include "GFace.h"
+#include "MElement.h"
+
+class elementFilter {
+public:
+  virtual bool operator () (MElement *) const = 0;
+};
+
+class elementFilterTrivial : public elementFilter {
+public:
+  bool operator () (MElement *) const {return true;}
+};
+
+class groupOfElements {
+public:
+  typedef std::set<MElement*> elementContainer; 
+  typedef std::set<MVertex*> vertexContainer; 
+private:
+  vertexContainer _vertices;
+  elementContainer _elements;
+  elementContainer _parents;
+public:
+  groupOfElements (int dim, int physical) {
+    addPhysical (dim, physical);
+  }   
+
+  groupOfElements (GFace*); 
+
+  void addPhysical(int dim, int physical) {
+    elementFilterTrivial filter;
+    addPhysical (dim, physical, filter);
+  }
+
+  void addElementary(GEntity *ge, const elementFilter &f);
+
+  void addPhysical(int dim, int physical, const elementFilter &);
+
+  vertexContainer::const_iterator vbegin () const {
+    return _vertices.begin();
+  }
+  vertexContainer::const_iterator vend () const {
+    return _vertices.end();
+  }
+  
+  elementContainer::const_iterator begin () const {
+    return _elements.begin();
+  }
+  elementContainer::const_iterator end () const {
+    return _elements.end();
+  }
+  size_t size () const {
+    return _elements.size();
+  }
+  // FIXME : NOT VERY ELEGANT !!!
+  bool find (MElement *e) const {
+    if (e->getParent() && _parents.find(e->getParent()) != _parents.end())return true;
+    return (_elements.find(e) != _elements.end()) ;
+  }
+  inline void insert (MElement *e) {
+    _elements.insert(e);
+    if (e->getParent()){
+      _parents.insert(e->getParent());      
+      for (int i=0;i<e->getParent()->getNumVertices();i++){
+	_vertices.insert(e->getParent()->getVertex(i));
+      }
+    }
+    else{
+      for (int i=0;i<e->getNumVertices();i++){
+	_vertices.insert(e->getVertex(i));
+      }
+    }
+  }
+};
+
+#endif
diff --git a/utils/api_demos/Makefile b/utils/api_demos/Makefile
index c6fc6009f74786834973cca2270daa1f9aebdcab..e05d66954b40b4dabb44c652207ed2a1fc2d3c50 100644
--- a/utils/api_demos/Makefile
+++ b/utils/api_demos/Makefile
@@ -4,6 +4,10 @@ GMSH = -L/usr/local/lib -lGmsh
 GLUT = -framework OpenGL -framework GLUT -framework Cocoa -framework ApplicationServices
 #GLUT = -lGLUT -lGLU -lGL -lX11
 
+mainElasticity: mainElasticity.cpp
+	g++ ${FLAGS} -o mainElasticity ${INC} mainElasticity.cpp\
+          ${GMSH} -llapack -lblas -lm -L/Users/remacle/SOURCES/gmsh/contrib/taucs/lib -ltaucs
+
 mainSimple: mainSimple.cpp
 	g++ ${FLAGS} -o mainSimple ${INC} mainSimple.cpp\
           ${GMSH} -llapack -lblas -lm
diff --git a/utils/api_demos/mainElasticity.cpp b/utils/api_demos/mainElasticity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..543a73b312c25f72341345e9d3334171356c7524
--- /dev/null
+++ b/utils/api_demos/mainElasticity.cpp
@@ -0,0 +1,32 @@
+#include <gmsh/Gmsh.h>
+#include <gmsh/elasticitySolver.h>
+#include <gmsh/PView.h>
+#include <gmsh/PViewData.h>
+
+int main (int argc, char* argv[]){
+  
+  if (argc != 2){
+    printf("usage : elasticity input_file_name\n");
+    return -1;
+  } 
+  
+  // globals are still present in Gmsh
+  GmshInitialize(argc, argv);
+  
+  // instanciate a solver
+  elasticitySolver mySolver (1000);
+  
+  // read some input file
+  mySolver.readInputFile(argv[1]);
+  
+  // solve the problem
+  mySolver.solve();
+
+  PView *pv = mySolver.buildDisplacementView("displacement");
+  pv->getData()->writeMSH("disp.msh", false);
+  delete pv;
+ 	
+  // stop gmsh
+  GmshFinalize();
+  
+}