diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index c67eb0e072416cdab3e2347fc8bd153b10507355..40bc19e3ffcbd5086813aa8c196a17dc2b252ae1 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -132,6 +132,28 @@ class dofManager{
   {
     numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
   }
+
+  inline void getDofValue(std::vector<Dof> &keys,std::vector<dataVec> &Vals)
+  {
+    int ndofs=keys.size();
+    Vals.reserve(Vals.size()+ndofs);
+    for (int i=0;i<ndofs;++i) Vals.push_back(getDofValue(keys[i]));
+  }
+
+  inline dataVec getDofValue(Dof key) const
+  {
+    {
+      typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
+      if (it != fixed.end()) return it->second;
+    }
+    {
+      std::map<Dof, int>::const_iterator it = unknown.find(key);
+      if (it != unknown.end())
+        return _current->getFromSolution(it->second);
+    }
+    return dataVec(0.0);
+  }
+
   inline dataVec getDofValue(int ent, int type) const
   {
     Dof key(ent, type);
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 2d2599c0b5575956f6302d6172fab4efcdedb261..e7e5575caee7d26bf71b5550f87648f53429e877 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -5,14 +5,12 @@
 
 #include <string.h>
 #include "GmshConfig.h"
-#include "elasticityTerm.h"
-#include "MTriangle.h"
-#include "MTetrahedron.h"
 #include "elasticitySolver.h"
 #include "linearSystemCSR.h"
 #include "linearSystemPETSc.h"
 #include "linearSystemGMM.h"
 #include "Numeric.h"
+#include "GModel.h"
 #include "functionSpace.h"
 #include "terms.h"
 #include "solverAlgorithms.h"
@@ -21,93 +19,17 @@
 #if defined(HAVE_POST)
 #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 (unsigned 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 (unsigned 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;  
-}
-
-#else
-PView* elasticitySolver::buildDisplacementView  (const std::string &postFileName)
-{
-  Msg::Error("Post-pro module not available");
-  return 0;
-}
 #endif
 
-
-static double vonMises(dofManager<double> *a, MElement *e, 
-                       double u, double v, double w, 
-                       double E, double nu, int _tag)
-{
-  double valx[256];
-  double valy[256];
-  double valz[256];
-  for (int k = 0; k < e->getNumVertices(); k++){
-    valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
-    valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
-    valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
-  }
-  double gradux[3];
-  double graduy[3];
-  double graduz[3];
-  e->interpolateGrad(valx, u, v, w, gradux);
-  e->interpolateGrad(valy, u, v, w, graduy);
-  e->interpolateGrad(valz, u, v, w, graduz);
-
-  double eps[6] = {gradux[0], graduy[1], graduz[2],
-                   0.5 * (gradux[1] + graduy[0]),
-                   0.5 * (gradux[2] + graduz[0]),
-                   0.5 * (graduy[2] + graduz[1])};
-  double A = E / (1. + nu);
-  double B = A * (nu / (1. - 2 * nu));
-  double trace = eps[0] + eps[1] + eps[2] ;
-  double sxx = A * eps[0] + B * trace;
-  double syy = A * eps[1] + B * trace;
-  double szz = A * eps[2] + B * trace;
-  double sxy = A * eps[3];
-  double sxz = A * eps[4];
-  double syz = A * eps[5];
-
-  double s[9] = {sxx, sxy, sxz,
-                 sxy, syy, syz,
-                 sxz, syz, szz};
-
-  return ComputeVonMises(s);
-}
-
 void elasticitySolver::setMesh(const std::string &meshFileName)
 {
   pModel = new GModel();
   pModel->readMSH(meshFileName.c_str());
   _dim = pModel->getNumRegions() ? 3 : 2;
+  if (LagSpace) delete LagSpace;
+  if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
+  if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+
 }
 
 void elasticitySolver::readInputFile(const std::string &fn)
@@ -116,13 +38,11 @@ void elasticitySolver::readInputFile(const std::string &fn)
   char what[256];
   while(!feof(f)){
     if(fscanf(f, "%s", what) != 1) return;
-    // printf("%s\n", what);
     if (!strcmp(what, "ElasticDomain")){
       elasticField field;
       int physical;
       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);
     }
@@ -133,75 +53,85 @@ void elasticitySolver::readInputFile(const std::string &fn)
       //      assign a tag
       //      elasticFields.push_back(field);
     }
-    else if (!strcmp(what, "NodalDisplacement")){
+    else if (!strcmp(what, "NodeDisplacement")){
       double val;
       int node, comp;
       if(fscanf(f, "%d %d %lf", &node, &comp, &val) != 3) return;
-      nodalDisplacements[ std::make_pair(node, comp) ] = val;
+      dirichletBC diri;
+      diri.g = new groupOfElements (0, node);
+      diri._f= simpleFunction<double>(val);
+      diri._comp=comp;
+      diri._tag=node;
+      diri.onWhat=BoundaryCondition::ON_VERTEX;
+      allDirichlet.push_back(diri);
     }
     else if (!strcmp(what, "EdgeDisplacement")){
       double val;
       int edge, comp;
       if(fscanf(f, "%d %d %lf", &edge, &comp, &val) != 3) return;
-      edgeDisplacements[ std::make_pair(edge, comp) ] = val;
       dirichletBC diri;
       diri.g = new groupOfElements (1, edge);
       diri._f= simpleFunction<double>(val);
       diri._comp=comp;
       diri._tag=edge;
-      edgeDiri.push_back(diri);
+      diri.onWhat=BoundaryCondition::ON_EDGE;
+      allDirichlet.push_back(diri);
     }
     else if (!strcmp(what, "FaceDisplacement")){
       double val;
       int face, comp;
       if(fscanf(f, "%d %d %lf", &face, &comp, &val) != 3) return;
-      faceDisplacements[ std::make_pair(face, comp) ] = val;
       dirichletBC diri;
       diri.g = new groupOfElements (2, face);
       diri._f= simpleFunction<double>(val);
       diri._comp=comp;
       diri._tag=face;
-      faceDiri.push_back(diri);    
+      diri.onWhat=BoundaryCondition::ON_FACE;
+      allDirichlet.push_back(diri);
     }
-    else if (!strcmp(what, "NodalForce")){
+    else if (!strcmp(what, "NodeForce")){
       double val1, val2, val3;
       int node;
       if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
-      nodalForces[node] = SVector3(val1, val2, val3);    
+      neumannBC neu;
+      neu.g = new groupOfElements (0, node);
+      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._tag=node;
+      neu.onWhat=BoundaryCondition::ON_VERTEX;
+      allNeumann.push_back(neu);  
     }
-    else if (!strcmp(what, "LineForce")){
+    else if (!strcmp(what, "EdgeForce")){
       double val1, val2, val3;
       int edge;
       if(fscanf(f, "%d %lf %lf %lf", &edge, &val1, &val2, &val3) != 4) return;
-      //printf("%d %lf %lf %lf\n", node, val1, val2, val3);
-      lineForces[edge] = SVector3(val1, val2, val3); 
       neumannBC neu;
       neu.g = new groupOfElements (1, edge);
       neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
       neu._tag=edge;
-      edgeNeu.push_back(neu);
+      neu.onWhat=BoundaryCondition::ON_EDGE;
+      allNeumann.push_back(neu);  
     }
     else if (!strcmp(what, "FaceForce")){
       double val1, val2, val3;
       int face;
       if(fscanf(f, "%d %lf %lf %lf", &face, &val1, &val2, &val3) != 4) return;
-      faceForces[face] = SVector3(val1, val2, val3);
       neumannBC neu;
       neu.g = new groupOfElements (2, face);
       neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
       neu._tag=face;
-      edgeNeu.push_back(neu);
+      neu.onWhat=BoundaryCondition::ON_FACE;
+      allNeumann.push_back(neu);  
     }
     else if (!strcmp(what, "VolumeForce")){
       double val1, val2, val3;
       int volume;
       if(fscanf(f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3) != 4) return;
-      volumeForces[volume] = SVector3(val1, val2, val3);
       neumannBC neu;
       neu.g = new groupOfElements (3, volume);
       neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
       neu._tag=volume;
-      edgeNeu.push_back(neu);  
+      neu.onWhat=BoundaryCondition::ON_VOLUME;
+      allNeumann.push_back(neu);  
     }
     else if (!strcmp(what, "MeshFile")){
       char name[245];
@@ -210,7 +140,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
     }
     else {
       Msg::Error("Invalid input : %s", what);
-      return;
+//      return;
     }
   }
   fclose(f);
@@ -226,236 +156,155 @@ void elasticitySolver::solve()
   linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
   lsys->setNoisy(2);
 #endif
+  if (pAssembler) delete pAssembler;
   pAssembler = new dofManager<double>(lsys);
 
-  std::map<int, std::vector<GEntity*> > groups[4];
-  pModel->getPhysicalGroups(groups);
-
   // we first do all fixations. the behavior of the dofManager is to 
   // give priority to fixations : when a dof is fixed, it cannot be
   // numbered afterwards
 
-  for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
-       it != nodalDisplacements.end(); ++it){
-    MVertex *v = pModel->getMeshVertexByTag(it->first.first);
-    pAssembler->fixVertex(v, it->first.second, _tag, it->second);
-    printf("-- Fixing node %3d comp %1d to %8.5f\n",
-           it->first.first, it->first.second, it->second);
-  }
-
-  for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
-       it != edgeDisplacements.end(); ++it){
-    elasticityTerm El(pModel, 1, 1, _tag);
-    El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag,
-                        simpleFunction<double>(it->second), *pAssembler);
-    printf("-- Fixing edge %3d comp %1d to %8.5f\n",
-           it->first.first, it->first.second, it->second);
-  }
 
-  for (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
-       it != faceDisplacements.end(); ++it){
-    elasticityTerm El(pModel, 1, 1, _tag);
-    El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag,
-                        simpleFunction<double>(it->second), *pAssembler);
-    printf("-- Fixing face %3d comp %1d to %8.5f\n",
-           it->first.first, it->first.second, it->second);
+  for (unsigned int i = 0; i < allDirichlet.size(); i++)
+  {
+    FilterDofComponent filter(allDirichlet[i]._comp);
+    if (allDirichlet[i].onWhat==BoundaryCondition::ON_VERTEX)
+      FixNodalDofs(*LagSpace,allDirichlet[i].g->vbegin(),allDirichlet[i].g->vend(),*pAssembler,allDirichlet[i]._f,filter);
+    else
+      FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
   }
 
-  // we number the nodes : when a node is numbered, it cannot be numbered
-  // again with another number. so we loop over elements
-  for (unsigned 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);	
-      }
-    }
+  // we number the dofs : when a dof is numbered, it cannot be numbered
+  // again with another number.
+  for (unsigned int i = 0; i < elasticFields.size(); ++i)
+  {
+    NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
   }
 
   // Now we start the assembly process
   // First build the force vector
-  // Nodes at geometrical vertices
-  for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
-       it != nodalForces.end(); ++it){
-    int iVertex = it->first;
-    SVector3 f = it->second;
-    std::vector<GEntity*> ent = groups[0][iVertex];
-    for (unsigned int i = 0; i < ent.size(); i++){
-      pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
-      pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
-      pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
-      printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
-             ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
-    }
-  }
 
-  // line forces
-  for (std::map<int, SVector3 >::iterator it = lineForces.begin();
-       it != lineForces.end(); ++it){
-    elasticityTerm El(pModel, 1, 1, _tag);
-    int iEdge = it->first;
-    SVector3 f = it->second;
-    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());
-  }
+  GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
 
-  // face forces
-  for (std::map<int, SVector3 >::iterator it = faceForces.begin();
-       it != faceForces.end(); ++it){
-    elasticityTerm El(pModel, 1, 1, _tag);
-    int iFace = it->first;
-    SVector3 f = it->second;
-    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());
+  for (unsigned int i = 0; i < allNeumann.size(); i++)
+  {
+    LoadTerm<FunctionSpace<SVector3> > Lterm(*LagSpace,allNeumann[i]._f);
+    if (allNeumann[i].onWhat==BoundaryCondition::ON_VERTEX)
+      Assemble(Lterm,*LagSpace,allNeumann[i].g->vbegin(),allNeumann[i].g->vend(),*pAssembler);
+    else
+      Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
   }
 
-  // volume forces
-  for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
-       it != volumeForces.end(); ++it){
-    elasticityTerm El(pModel, 1, 1, _tag);
-    int iVolume = it->first;
-    SVector3 f = it->second;
-    El.setVector(f);
-    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++){
-      //   to do
-      //      El.addToRightHandSide(*pAssembler, ent[i]);
-    }
-  }
+// bulk material law
 
-  // assembling the stifness matrix
-  for (unsigned int i = 0; i < elasticFields.size(); i++){
-    SElement::setShapeEnrichement(elasticFields[i]._enrichment);
-    for (unsigned 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);      
-    }
+  GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+  for (unsigned int i = 0; i < elasticFields.size(); i++)
+  {
+    IsotropicElasticTerm<FunctionSpace<SVector3> ,FunctionSpace<SVector3> > Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
+    Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
   }
-/*
+
   printf("-- done assembling!\n");
-    for (int i=0;i<pAssembler->sizeOfR();i++)
-    {
-        if (fabs(lsys->getFromRightHandSide(i))>0.000001)
-        printf("%g ",lsys->getFromRightHandSide(i));
-    }
-    printf("\n");
-*/
-  // solving
   lsys->systemSolve();
   printf("-- done solving!\n");
 }
 
 
 
+#if defined(HAVE_POST)
 
-void MyelasticitySolver::solve()
+static double vonMises(dofManager<double> *a, MElement *e, 
+                       double u, double v, double w, 
+                       double E, double nu, int _tag)
 {
-#if defined(HAVE_TAUCS)
-  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
-#elif defined(HAVE_PETSC)
-  linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
-#else
-  linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
-  lsys->setNoisy(2);
-#endif
-  pAssembler = new dofManager<double>(lsys);
-
-  std::map<int, std::vector<GEntity*> > groups[4];
-  pModel->getPhysicalGroups(groups);
-
-  // we first do all fixations. the behavior of the dofManager is to 
-  // give priority to fixations : when a dof is fixed, it cannot be
-  // numbered afterwards
-  VectorLagrangeFunctionSpace P123(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
-
-  for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
-       it != nodalDisplacements.end(); ++it){
-    MVertex *v = pModel->getMeshVertexByTag(it->first.first);
-    pAssembler->fixVertex(v, it->first.second, _tag, it->second);
-    printf("-- Fixing node %3d comp %1d to %8.5f\n",
-           it->first.first, it->first.second, it->second);
+  double valx[256];
+  double valy[256];
+  double valz[256];
+  for (int k = 0; k < e->getNumVertices(); k++){
+    valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
+    valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
+    valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
   }
+  double gradux[3];
+  double graduy[3];
+  double graduz[3];
+  e->interpolateGrad(valx, u, v, w, gradux);
+  e->interpolateGrad(valy, u, v, w, graduy);
+  e->interpolateGrad(valz, u, v, w, graduz);
 
+  double eps[6] = {gradux[0], graduy[1], graduz[2],
+                   0.5 * (gradux[1] + graduy[0]),
+                   0.5 * (gradux[2] + graduz[0]),
+                   0.5 * (graduy[2] + graduz[1])};
+  double A = E / (1. + nu);
+  double B = A * (nu / (1. - 2 * nu));
+  double trace = eps[0] + eps[1] + eps[2] ;
+  double sxx = A * eps[0] + B * trace;
+  double syy = A * eps[1] + B * trace;
+  double szz = A * eps[2] + B * trace;
+  double sxy = A * eps[3];
+  double sxz = A * eps[4];
+  double syz = A * eps[5];
 
-  for (unsigned int i = 0; i < edgeDiri.size(); i++)
-  {
-    FixNodalDofs(P123,edgeDiri[i].g->begin(),edgeDiri[i].g->end(),*pAssembler,edgeDiri[i]._f,FilterDofComponent(edgeDiri[i]._comp));
-    printf("-- Fixing edge %3d comp %3d \n", edgeDiri[i]._tag, edgeDiri[i]._comp);
-  }
+  double s[9] = {sxx, sxy, sxz,
+                 sxy, syy, syz,
+                 sxz, syz, szz};
 
-  for (unsigned int i = 0; i < faceDiri.size(); i++)
-  {
-    FixNodalDofs(P123,faceDiri[i].g->begin(),faceDiri[i].g->end(),*pAssembler,faceDiri[i]._f,FilterDofComponent(faceDiri[i]._comp));
-    printf("-- Fixing face %3d comp %3d \n", faceDiri[i]._tag, faceDiri[i]._comp);
-  }
+  return ComputeVonMises(s);
+}
 
-  // we number the nodes : when a node is numbered, it cannot be numbered
-  // again with another number. so we loop over elements
+PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
+{
+  std::set<MVertex*> v;
   for (unsigned int i = 0; i < elasticFields.size(); ++i)
   {
-    NumberDofs(P123, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
-  }
-
-  // Now we start the assembly process
-  // First build the force vector
-  // Nodes at geometrical vertices
-  for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
-       it != nodalForces.end(); ++it){
-    int iVertex = it->first;
-    SVector3 f = it->second;
-    std::vector<GEntity*> ent = groups[0][iVertex];
-    for (unsigned int i = 0; i < ent.size(); i++){
-      pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
-      pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
-      pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
-      printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
-             ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
+    for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
+    {
+      MElement *e=*it;
+      for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
     }
   }
-
-
-  GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
-
-  for (unsigned int i = 0; i < edgeNeu.size(); i++)
+  std::map<int, std::vector<double> > data;  
+  for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
   {
-    LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,edgeNeu[i]._f);
-    Assemble(Lterm,P123,edgeNeu[i].g->begin(),edgeNeu[i].g->end(),Integ_Boundary,*pAssembler);
-    printf("-- Force on edge %3d \n", edgeNeu[i]._tag);
+    SVector3 val(0,0,0);
+    for (unsigned int i = 0; i < elasticFields.size(); ++i)
+    {
+      std::vector<Dof> D;
+      std::vector<SVector3> SFVals;
+      std::vector<double> Vals;
+      LagSpace->getKeys(*it,D);
+      pAssembler->getDofValue(D,Vals);
+      LagSpace->f(*it,SFVals);
+      for (int i=0;i<D.size();++i)
+        val+=SFVals[i]*Vals[i];
+    }
+    std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
+    data[(*it)->getNum()]=vec;
   }
+  PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
+  return pv;  
+}
 
-  for (unsigned int i = 0; i < faceNeu.size(); i++)
-  {
-    LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,faceNeu[i]._f);
-    Assemble(Lterm,P123,faceNeu[i].g->begin(),faceNeu[i].g->end(),Integ_Boundary,*pAssembler);
-    printf("-- Force on face %3d \n", faceNeu[i]._tag);
-  }
+PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
+{
+  std::map<int, std::vector<double> > data;  
+  PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
+  return pv;  
+}
 
-  for (unsigned int i = 0; i < volumeNeu.size(); i++)
-  {
-    LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,volumeNeu[i]._f);
-    Assemble(Lterm,P123,volumeNeu[i].g->begin(),volumeNeu[i].g->end(),Integ_Boundary,*pAssembler);
-    printf("-- Force on volume %3d \n", volumeNeu[i]._tag);
-  }
 
-  GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
-  for (unsigned int i = 0; i < elasticFields.size(); i++)
-  {
-    IsotropicElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
-    Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
-  }
+#else
+PView* elasticitySolver::buildDisplacementView  (const std::string &postFileName)
+{
+  Msg::Error("Post-pro module not available");
+  return 0;
+}
 
-  printf("-- done assembling!\n");
-  lsys->systemSolve();
-  printf("-- done solving!\n");
+PView* elasticitySolver::buildVonMisesView(const std::string &postFileName)
+{
+  Msg::Error("Post-pro module not available");
+  return 0;
 }
+
+#endif
+
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 5a97d9cbb45adba80302732b4d193b399db72ab9..d9b82a8fe77b09495bf6c1c4a46c5949601e297c 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -11,6 +11,7 @@
 #include "SVector3.h"
 #include "dofManager.h"
 #include "simpleFunction.h"
+#include "functionSpace.h"
 
 class GModel;
 class PView;
@@ -24,68 +25,56 @@ struct elasticField {
   elasticField () : g(0), _enrichment(0),_tag(0){}
 };
 
-struct dirichletBC {
+struct BoundaryCondition
+{
+  enum location{UNDEF,ON_VERTEX,ON_EDGE,ON_FACE,ON_VOLUME};
+  location onWhat; // on vertices or elements
   int _tag; // tag for the dofManager
-  int _comp; // component
   groupOfElements *g; // support for this BC
+  BoundaryCondition() : g(0),_tag(0),onWhat(UNDEF) {};
+};
+
+struct dirichletBC : public BoundaryCondition
+{
+  int _comp; // component
   simpleFunction<double> _f;
-  dirichletBC () : g(0),_comp(0),_tag(0){}
+  dirichletBC () :BoundaryCondition(),_comp(0),_f(0){}
 };
 
-struct neumannBC {
-  int _tag; // tag for the dofManager
-  groupOfElements *g; // support for this BC
+struct neumannBC  : public BoundaryCondition
+{
   simpleFunction<SVector3> _f;
-  neumannBC () : g(0),_tag(0),_f(SVector3(0,0,0)){}
+  neumannBC () : BoundaryCondition(),_f(SVector3(0,0,0)){}
 };
 
 // an elastic solver ...
-class elasticitySolver{
+class elasticitySolver
+{
  protected:
   GModel *pModel;
   int _dim, _tag;
   dofManager<double> *pAssembler;
+  FunctionSpace<SVector3> *LagSpace;
+
   // young modulus and poisson coefficient per physical
   std::vector<elasticField> elasticFields;
-  // imposed nodal forces
-  std::map<int, SVector3> nodalForces;
-  // imposed line forces
-  std::map<int, SVector3> lineForces;
-  std::vector<neumannBC> edgeNeu;
-  // imposed face forces
-  std::map<int, SVector3> faceForces;
-  std::vector<neumannBC> faceNeu;
-  // imposed volume forces
-  std::map<int, SVector3> volumeForces;
-  std::vector<neumannBC> volumeNeu;
-  // imposed nodal displacements
-  std::map<std::pair<int,int>, double> nodalDisplacements;
-  // imposed edge displacements
-  std::map<std::pair<int,int>, double> edgeDisplacements;
-  std::vector<dirichletBC> edgeDiri;
-  // imposed face displacements
-  std::map<std::pair<int,int>, double> faceDisplacements;
-  std::vector<dirichletBC> faceDiri;
-//  std::vector<dirichletBC> faceDiri;
+  // neumann BC
+  std::vector<neumannBC> allNeumann;
+  // dirichlet BC
+  std::vector<dirichletBC> allDirichlet;
+  
  public:
-  elasticitySolver(int tag) : _tag(tag) {}
-  void addNodalForces (int iNode, const SVector3 &f)
+  elasticitySolver(int tag) : _tag(tag),LagSpace(0),pAssembler(0) {}
+  virtual ~elasticitySolver()
   {
-    nodalForces[iNode] = f;
+    if (LagSpace) delete LagSpace;
+    if (pAssembler) delete pAssembler;
   }
-  void addNodalDisplacement(int iNode, int dir, double val)
-  {
-    nodalDisplacements[std::make_pair(iNode, dir)] = val;
-  }
-  //  void addElasticConstants(double e, double nu, int physical)
-  //  {
-  //    elasticConstants[physical] = std::make_pair(e, nu);
-  //  }
+  void readInputFile(const std::string &meshFileName);
   void setMesh(const std::string &meshFileName);
   virtual void solve();
-  void readInputFile(const std::string &meshFileName);
   virtual PView *buildDisplacementView(const std::string &postFileName);
-  // PView *buildVonMisesView(const std::string &postFileName);
+  virtual PView *buildVonMisesView(const std::string &postFileName);
   // std::pair<PView *, PView*> buildErrorEstimateView
   //   (const std::string &errorFileName, double, int);
   // std::pair<PView *, PView*> buildErrorEstimateView
@@ -93,15 +82,4 @@ class elasticitySolver{
 };
 
 
-
-// another elastic solver ...
-class MyelasticitySolver : public elasticitySolver
-{
- public:
-  MyelasticitySolver(int tag) : elasticitySolver(tag) {}
-  virtual void solve();
-};
-
-
-
 #endif
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 105fea2127d7129472567354f7cadcad2e1e8457..659037c2f4a3414d72a0648c25d704b8e3658de8 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -67,13 +67,15 @@ class FunctionSpaceBase
 {
  public:
   virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
-  virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0; 
+  virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+  virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
+  virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
 };
 
 template<class T>
 class FunctionSpace : public FunctionSpaceBase
 {
- protected:
+ public:
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
 /*  typedef typename TensorialTraits<T>::HessType HessType;
@@ -81,6 +83,7 @@ class FunctionSpace : public FunctionSpaceBase
   typedef typename TensorialTraits<T>::CurlType CurlType;*/
  
  public:
+  virtual int f(MVertex *ver, std::vector<ValType> &vals) {} // used for neumann BC
   virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0;
   virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)=0;
 //  virtual groupOfElements* getSupport()=0;// probablement inutile
@@ -89,26 +92,24 @@ class FunctionSpace : public FunctionSpaceBase
   virtual int divf(MElement *ele, double u, double v, double w,std::vector<DivType> &divs);
   virtual int curlf(MElement *ele, double u, double v, double w,std::vector<CurlType> &curls);*/
   virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
-  virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0; 
+  virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+  virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
+  virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
 };
 
 class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
 {
- protected:
+ public:
   typedef TensorialTraits<double>::ValType ValType;
   typedef TensorialTraits<double>::GradType GradType;
   typedef TensorialTraits<double>::HessType HessType;
   typedef TensorialTraits<double>::DivType DivType;
   typedef TensorialTraits<double>::CurlType CurlType;
+ protected:
   std::vector<basisFunction*> basisFunctions;
 
   int _iField; // field number (used to build dof keys)
 
-  Dof getLocalDof(MElement *ele, int i) const
-  {
-    return Dof(ele->getVertex(i)->getNum(), _iField);
-  }
-
  public:
   ScalarLagrangeFunctionSpace(int i=0):_iField(i) {}
   virtual int getId(void) const {return _iField;};
@@ -118,7 +119,8 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
     int curpos=vals.size();
     vals.resize(curpos+ndofs);
     ele->getShapeFunctions(u, v, w, &(vals[curpos]));
-  }; 
+  };
+  virtual int f(MVertex *ver, std::vector<ValType> &vals) {vals.push_back(1.0);} // used for neumann BC
   virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
   {
     int ndofs= ele->getNumVertices();
@@ -143,8 +145,13 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
       int ndofs= ele->getNumVertices();
       keys.reserve(keys.size()+ndofs);
       for (int i=0;i<ndofs;++i)
-        keys.push_back(getLocalDof(ele,i));
+        getKeys(ele->getVertex(i),keys);
   }
+  virtual int getNumKeys(MVertex *ver) { return 1;}
+  virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+  {
+    keys.push_back(Dof(ver->getNum(), _iField));
+  };
 };
 
 
@@ -180,6 +187,20 @@ public :
   }
 
   virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;}
+  virtual int f(MVertex *ver, std::vector<ValType> &vals)
+  {
+    std::vector<double> valsd;
+    ScalarFS->f(ver,valsd);
+    int nbdofs=valsd.size();
+    int nbcomp=comp.size();
+    int curpos=vals.size();
+    vals.reserve(curpos+nbcomp*nbdofs);
+    for (int j=0;j<nbcomp;++j)
+    {
+      for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]);
+    }
+  } // used for neumann BC
+
   virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
   {
     std::vector<double> valsd;
@@ -236,14 +257,36 @@ public :
       }
     }
   }
+
+  virtual int getNumKeys(MVertex *ver) { return ScalarFS->getNumKeys(ver)*comp.size();}
+  virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+  {
+    int nk=ScalarFS->getNumKeys(ver);
+    std::vector<Dof> bufk;
+    bufk.reserve(nk);
+    ScalarFS->getKeys(ver,bufk);
+    int nbdofs=bufk.size();
+    int nbcomp=comp.size();
+    int curpos=keys.size();
+    keys.reserve(curpos+nbcomp*nbdofs);
+    for (int j=0;j<nbcomp;++j)
+    {
+      for (int i=0;i<nk;++i)
+      {
+        int i1,i2;
+        Dof::getTwoIntsFromType(bufk[i].getType(), i1,i2);
+        keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i1)));
+      }
+    }
+  };
 };
 
 class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
 {
+ protected:
+  static const SVector3 BasisVectors[3];
  public:
   enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 };
-  static const SVector3 BasisVectors[3];
-
   typedef TensorialTraits<SVector3>::ValType ValType;
   typedef TensorialTraits<SVector3>::GradType GradType;
   typedef TensorialTraits<SVector3>::HessType HessType;
@@ -305,6 +348,12 @@ class CompositeFunctionSpace : public FunctionSpace<T>
       delete (*it);
   }
 
+  virtual int f(MVertex *ver, std::vector<ValType> &vals)
+  {
+    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+      (*it)->f(ver,vals);
+  }
+
   virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
   {
     for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
@@ -330,11 +379,26 @@ class CompositeFunctionSpace : public FunctionSpace<T>
     for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
       (*it)->getKeys(ele,keys);
   }
+  
+  virtual int getNumKeys(MVertex *ver)
+  {
+    int ndofs=0;
+    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+      ndofs+=(*it)->getNumKeys(ver);
+    return ndofs;
+  }
+
+  virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+  {
+    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+      (*it)->getKeys(ver,keys);
+  };
 };
 
 template<class T>
 class xFemFunctionSpace : public FunctionSpace<T>
 {
+ public:
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
   typedef typename TensorialTraits<T>::HessType HessType;
@@ -344,10 +408,14 @@ class xFemFunctionSpace : public FunctionSpace<T>
   FunctionSpace<T>* _spacebase;
 //  Function<double>* enrichment;
  public:
+  virtual int f(MVertex *ver, std::vector<ValType> &vals);
   virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
   virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
-  int getNumKeys(MElement *ele);
-  void getKeys(MElement *ele, std::vector<Dof> &keys);
+  virtual int getNumKeys(MElement *ele);
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
+  virtual int getNumKeys(MVertex *ver);
+  virtual int getKeys(MVertex *ver, std::vector<Dof> &keys);
+
 };
 
 
@@ -364,10 +432,13 @@ class FilteredFunctionSpace : public FunctionSpace<T>
   F &_filter;
   
  public:
+  virtual int f(MVertex *ver, std::vector<ValType> &vals);
   virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
   virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
-  int getNumKeys(MElement *ele);
-  void getKeys(MElement *ele, std::vector<Dof> &keys);
+  virtual int getNumKeys(MElement *ele);
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
+  virtual int getNumKeys(MVertex *ver);
+  virtual int getKeys(MVertex *ver, std::vector<Dof> &keys);
 };
 
 
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 70c1abb84860550eb746eeb4fe194641aad9efb5..0afc37ef5120072056f8e6c20fe5b3e74174c83d 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -21,6 +21,7 @@
 #include "MVertex.h"
 
 
+
 template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler) // symmetric
 {
   fullMatrix<typename Assembler::dataMat> localMatrix;
@@ -64,6 +65,20 @@ template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,Func
   }
 }
 
+template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
+{
+  fullVector<typename Assembler::dataMat> localVector;
+  std::vector<Dof> R;
+  for (Iterator it = itbegin;it!=itend; ++it)
+  {
+    MVertex *v = *it;
+    R.clear();
+    term.get(v,localVector);
+    space.getKeys(v,R);
+    assembler.assemble(R, localVector);
+  }
+}
+
 template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler)
 {
   fullVector<typename Assembler::dataMat> localVector;
@@ -75,7 +90,6 @@ template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &
   assembler.assemble(R, localVector);
 }
 
-
 template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &dofs,std::vector<typename Assembler::dataVec> &vals)
 {
   int nbff=dofs.size();
@@ -86,6 +100,12 @@ template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &do
 }
 
 class FilterDof
+{
+ public:
+  virtual bool operator()(Dof key)=0;
+};
+
+class FilterDofTrivial
 {
  public:
   virtual bool operator()(Dof key) {return true;}
@@ -100,43 +120,60 @@ class FilterDofComponent :public FilterDof
   {
     int type=key.getType();
     int icomp,iphys;
-    Dof::getTwoIntsFromType(type, iphys, icomp);
+    Dof::getTwoIntsFromType(type, icomp, iphys);
     if (icomp==comp) return true;
     return false;
   }
 };
 
 
-template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof filter)
+template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
 {
-  for (Iterator it=itbegin;it!=itend;++it)
+  std::vector<MVertex*> tabV;
+  int nv=e->getNumVertices();
+  std::vector<Dof> R;
+  space.getKeys(e,R);
+  tabV.reserve(nv);
+  for (int i=0;i<nv;++i) tabV.push_back(e->getVertex(i));
+  
+  for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
   {
-    MElement *e=*it;
-    std::vector<MVertex*> tabV;
-    int nv=e->getNumVertices();
-    std::vector<Dof> R;
-    space.getKeys(e,R);
-    tabV.reserve(nv);
-    for (int i=0;i<nv;++i) tabV.push_back(e->getVertex(i));
-    
-    for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
+    Dof key=*itd;
+    if (filter(key))
     {
-      Dof key=*itd;
-      if (filter(key))
+      for (int i=0;i<nv;++i)
       {
-        for (int i=0;i<nv;++i)
+        if (tabV[i]->getNum()==key.getEntity())
         {
-          if (tabV[i]->getNum()==key.getEntity())
-          {
-            assembler.fixDof(key, fct(tabV[i]->x(),tabV[i]->y(),tabV[i]->z()));
-            break;
-          }
+          assembler.fixDof(key, fct(tabV[i]->x(),tabV[i]->y(),tabV[i]->z()));
+          break;
         }
       }
     }
   }
 }
 
+template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MVertex *v,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
+{
+  std::vector<Dof> R;
+  space.getKeys(v,R);
+  for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
+  {
+    Dof key=*itd;
+    if (filter(key))
+    {
+      assembler.fixDof(key, fct(v->x(),v->y(),v->z()));
+    }
+  }
+}
+
+template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
+{
+  for (Iterator it=itbegin;it!=itend;++it)
+  {
+    FixNodalDofs(space,*it,assembler,fct,filter);
+  }
+}
 
 template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
 {
diff --git a/Solver/terms.h b/Solver/terms.h
index 721574d943ab6c014b66654e3129503a51cce76a..b5734be749fba741b64dbafdb670a3e57bf96e0c 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -22,6 +22,7 @@
 #include "groupOfElements.h"
 
 
+
 // evaluation of a field ???
 template<class T>
 class Field {
@@ -80,6 +81,7 @@ class  LinearTermBase
 {
   public:
   virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
+  virtual void get(MVertex *v,fullVector<double> &m) =0;
 };
 
 template<class S1> class LinearTerm : public LinearTermBase
@@ -94,7 +96,7 @@ template<class S1> class LinearTerm : public LinearTermBase
 class  ScalarTermBase
 {
  public : 
-  virtual void get(MElement *ele,int npts,IntPt *GP,double &v) =0;
+  virtual void get(MElement *ele,int npts,IntPt *GP,double &val) =0;
 };
 
 class ScalarTerm : public ScalarTermBase
@@ -105,6 +107,60 @@ class ScalarTerm : public ScalarTermBase
 
 
 
+template<class S1,class S2> class LaplaceTerm : public BilinearTerm<S1,S2> 
+{
+ public : 
+  LaplaceTerm(S1& space1_,S2& space2_) : BilinearTerm<S1,S2>(space1_,space2_)
+  {}
+  
+  virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+  {
+    Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+  }
+  virtual void get(MVertex *v,fullMatrix<double> &m)
+  {
+    Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+  }
+}; // class
+
+
+
+template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
+{
+ public : 
+  LaplaceTerm(S1& space1_) : BilinearTerm<S1,S1>(space1_,space1_)
+  {}
+  
+  virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+  {
+    int nbFF = BilinearTerm<S1,S1>::space1.getNumKeys(ele);
+    double jac[3][3];
+    m.resize(nbFF, nbFF);
+    m.setAll(0.);
+    for (int i = 0; i < npts; i++)
+    {
+      const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
+      const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+      std::vector<typename S1::GradType> Grads;
+      BilinearTerm<S1,S1>::space1.gradf(ele,u, v, w, Grads);
+      for (int j = 0; j < nbFF; j++)
+      {
+        for (int k = j; k < nbFF; k++)
+        {
+          double contrib=weight * detJ * dot(Grads[j],Grads[k]);
+          m(j,k)+=contrib;
+          if (j!=k) m(k,j)+=contrib;
+        }
+      }
+    }
+//    m.print("");
+//    exit(0);
+  }
+}; // class
+
+
+
+
 template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,S2> // non symmetric 
 {
  protected : 
@@ -238,9 +294,9 @@ inline double dot(const double &a, const double &b)
 
 template<class S1> class LoadTerm : public LinearTerm<S1>
 {
-  simpleFunction<typename S1::ValType> Load;
+  simpleFunction<typename S1::ValType> &Load;
  public : 
-  LoadTerm(S1& space1_,simpleFunction<typename S1::ValType> Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
+  LoadTerm(S1& space1_,simpleFunction<typename S1::ValType> &Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
   virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
   {
     double nbFF=LinearTerm<S1>::space1.getNumKeys(ele);
@@ -253,12 +309,29 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
       const double weight = GP[i].weight;const double detJ = ele->getJacobian(u, v, w, jac);
       std::vector<typename S1::ValType> Vals;
       LinearTerm<S1>::space1.f(ele,u, v, w, Vals);
+      SPoint3 p;
+      ele->pnt(u, v, w, p);
+      typename S1::ValType load=Load(p.x(),p.y(),p.z());
       for (int j = 0; j < nbFF ; ++j)
       {
-        m(j)+=dot(Vals[j],Load(u,v,w))*weight*detJ;
+        m(j)+=dot(Vals[j],load)*weight*detJ;
       }
     }
   }
+
+  virtual void get(MVertex *ver,fullVector<double> &m)
+  {
+    double nbFF=LinearTerm<S1>::space1.getNumKeys(ver);
+    double jac[3][3];
+    m.resize(nbFF);
+    std::vector<typename S1::ValType> Vals;
+    LinearTerm<S1>::space1.f(ver, Vals);
+    typename S1::ValType load=Load(ver->x(),ver->y(),ver->z());
+    for (int j = 0; j < nbFF ; ++j)
+    {
+      m(j)=dot(Vals[j],load);
+    }
+  }
 };