diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 1b7880f3f4960de7960b8aaad472abc32c630dfd..d45a69f6076756bca7c3e4291c1864c351a214d1 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -53,7 +53,8 @@ class FunctionSpace : public FunctionSpaceBase
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
   virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0;
-  virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)=0;
+  virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) =0;
+  virtual void gradfuvw(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) {}; // should return to pure virtual once all is done.
   virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0;
 };
@@ -102,6 +103,20 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
                             ));
   }
 
+  virtual void gradfuvw(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+  {
+    int ndofs= ele->getNumVertices();
+    grads.reserve(grads.size()+ndofs);
+    double gradsuvw[256][3];
+    ele->getGradShapeFunctions(u, v, w, gradsuvw);
+    for (int i=0;i<ndofs;++i)
+      grads.push_back(GradType(
+      gradsuvw[i][0] ,
+      gradsuvw[i][1] ,
+      gradsuvw[i][2]
+                      ));
+  }
+
   virtual int getNumKeys(MElement *ele) {return ele->getNumVertices();}
 
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
@@ -177,6 +192,26 @@ public :
     }
   }
   
+  
+  virtual void gradfuvw(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+  {
+    std::vector<SVector3> gradsd;
+    ScalarFS->gradfuvw(ele,u,v,w,gradsd);
+    int nbdofs=gradsd.size();
+    int nbcomp=comp.size();
+    int curpos=grads.size();
+    grads.reserve(curpos+nbcomp*nbdofs);
+    GradType val;
+    for (int j=0;j<nbcomp;++j)
+    {
+      for (int i=0;i<nbdofs;++i)
+      {
+        tensprod(multipliers[j],gradsd[i],val);
+        grads.push_back(val);
+      }
+    }
+  }
+
   virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele)*comp.size();}
   
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys)
diff --git a/Solver/groupOfElements.cpp b/Solver/groupOfElements.cpp
index 74cd63e139bc4b6e24ee9e30618ae8b2920ae738..a0984909cb836a5e67fcbcc1fbb21a60b858c00a 100644
--- a/Solver/groupOfElements.cpp
+++ b/Solver/groupOfElements.cpp
@@ -21,7 +21,7 @@ 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];
+  std::vector<GEntity*> &ent = groups[dim][physical];
   for (unsigned int i = 0; i < ent.size(); i++){
     addElementary(ent[i], filter);
   }