diff --git a/NonLinearSolver/Domain/partDomain.cpp b/NonLinearSolver/Domain/partDomain.cpp
index a02b0b9bbcd275ea413d445d70b26b95290ce54e..b5d4483ae0d07644e063f221f9917891c2866c3d 100644
--- a/NonLinearSolver/Domain/partDomain.cpp
+++ b/NonLinearSolver/Domain/partDomain.cpp
@@ -428,7 +428,7 @@ void partDomain::computeAverageHighOrderStress(const IPField* ipf, STensor33& si
 	if (g->size()>0){
 		IntPt* GP;
 		const ipFiniteStrain *vipv[256];
-		FunctionSpaceBase* space = this->getFunctionSpace();
+		const FunctionSpaceBase* space = this->getFunctionSpace();
 		for (groupOfElements::elementContainer::iterator it= g->begin(); it!=g->end(); it++){
 			MElement* ele= *it;
 			if(ele->getParent()) ele = ele->getParent();
diff --git a/NonLinearSolver/Domain/partDomain.h b/NonLinearSolver/Domain/partDomain.h
index 602f82ecc5b3d61bcf12009e7bef734b78d37e69..9c41edc39216d466aa48ff72e63a68ec44df36f6 100644
--- a/NonLinearSolver/Domain/partDomain.h
+++ b/NonLinearSolver/Domain/partDomain.h
@@ -188,8 +188,9 @@ public:
   // true is domain has interface terms
   virtual bool IsInterfaceTerms() const=0;
   functionSpaceType::whichSpace getFunctionSpaceType() const{return _wsp;}
-  // can be return const FunctionSpaceBase if the function of this class are declarated const
-  virtual FunctionSpaceBase* getFunctionSpace() const=0;
+
+  virtual FunctionSpaceBase* getFunctionSpace() =0;
+  virtual const FunctionSpaceBase* getFunctionSpace() const=0;
 // some data of BC have to be set by domain
   virtual FunctionSpaceBase* getSpaceForBC(const nonLinearBoundaryCondition::type bc_type, const nonLinearBoundaryCondition::location bc_location,
                                            const int dof_comp,const groupOfElements *groupBC) const=0; // for dirichlet, neumann and initial
@@ -415,10 +416,18 @@ class dgPartDomain : public partDomain{
   virtual void createIPState(AllIPState::ipstateContainer& map, const bool *state) const;
   virtual int getDim() const=0;
   virtual bool IsInterfaceTerms() const{return true;}
-  virtual FunctionSpaceBase* getFunctionSpace() const=0;
-  virtual FunctionSpaceBase* getFunctionSpaceMinus() const=0;
-  virtual FunctionSpaceBase* getFunctionSpacePlus() const=0;
-  virtual FunctionSpaceBase* getInterfaceFunctionSpace() const =0;
+  virtual FunctionSpaceBase* getFunctionSpace() =0;
+  virtual const FunctionSpaceBase* getFunctionSpace() const=0;
+  
+  virtual FunctionSpaceBase* getFunctionSpaceMinus()=0;
+  virtual const FunctionSpaceBase* getFunctionSpaceMinus() const=0;
+  
+  virtual FunctionSpaceBase* getFunctionSpacePlus()=0;
+  virtual const FunctionSpaceBase* getFunctionSpacePlus() const=0;
+  
+  virtual FunctionSpaceBase* getInterfaceFunctionSpace() =0;
+  virtual const FunctionSpaceBase* getInterfaceFunctionSpace() const =0;
+  
   virtual void matrixByPerturbation(const int ibulk, const int iinter, const int ivirt,const double eps=1e-8)=0;
   virtual void setMaterialLaw(const std::map<int,materialLaw*> &maplaw)=0;
   virtual materialLaw* getMaterialLaw(){Msg::Error("The law to retrieve is not given on a dgdom"); return NULL;}
diff --git a/NonLinearSolver/contact/contactFunctionSpace.h b/NonLinearSolver/contact/contactFunctionSpace.h
index 19a8e31d285d0285ba9b855c8742184c6a3c8a4e..fc775177b69687fda0f6085d36daaa69c8ef5ed2 100644
--- a/NonLinearSolver/contact/contactFunctionSpace.h
+++ b/NonLinearSolver/contact/contactFunctionSpace.h
@@ -34,10 +34,10 @@ template <class T1> class ContactSpaceBase : public FunctionSpaceBase{
   }
 	virtual int getId(void) const{return _id;};
 
-  virtual void getKeys(MElement *ele,std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele,std::vector<Dof> &keys) const{
     _space->getKeys(ele,keys);
   }
-  virtual int getNumKeys(MElement *ele){
+  virtual int getNumKeys(MElement *ele) const{
     return _space->getNumKeys(ele);
   }
 };
@@ -49,12 +49,12 @@ class rigidContactSpaceBase : public ContactSpaceBase<double>{
   rigidContactSpaceBase(const int id, FunctionSpace<double> *sp, MVertex *ver) : ContactSpaceBase<double>(id,sp), _gc(ver){}
   rigidContactSpaceBase(const int id, FunctionSpace<double> *sp, MVertex *ver, const int comp1,
                    const int comp2 = -1, const int comp3 =-1) : ContactSpaceBase<double>(id,sp,comp1,comp2,comp3),_gc(ver){}
-  virtual void getKeys(MElement *ele,std::vector<Dof> &keys)=0;
-  virtual int getNumKeys(MElement *ele)=0;
-  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys)=0;
-  virtual int getNumKeysOfGravityCenter()=0;
-  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys)=0;
-  virtual int getNumKeys(MElement *ele, int ind)=0; // number key in one component + number of key of GC
+  virtual void getKeys(MElement *ele,std::vector<Dof> &keys) const=0;
+  virtual int getNumKeys(MElement *ele) const =0;
+  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys) const=0;
+  virtual int getNumKeysOfGravityCenter() const=0;
+  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys) const=0;
+  virtual int getNumKeys(MElement *ele, int ind) const=0; // number key in one component + number of key of GC
 };
 #endif // CONTACTSPACE_H_
 
diff --git a/NonLinearSolver/field/energyField.cpp b/NonLinearSolver/field/energyField.cpp
index 5c43836e0b38e37881f3cd24dbf182848a3f42bc..f307d8f861a58be6c41a8ff118aa6d344ec78e23 100644
--- a/NonLinearSolver/field/energyField.cpp
+++ b/NonLinearSolver/field/energyField.cpp
@@ -175,7 +175,7 @@ double energeticField::kineticEnergy(MElement *ele, const partDomain *dom) const
     std::vector<Dof> R;
     std::vector<double> velocities;
     std::vector<double> vmass;
-    FunctionSpaceBase *sp = dom->getFunctionSpace();
+    const FunctionSpaceBase *sp = dom->getFunctionSpace();
     sp->getKeys(ele,R);
     int nkeys = sp->getNumKeys(ele);
     _solver->getStaticDofManager()->getDofValue(R,velocities,nonLinearBoundaryCondition::velocity);
diff --git a/NonLinearSolver/internalPoints/IntPtData.h b/NonLinearSolver/internalPoints/IntPtData.h
index ca4f4e520556567a424e2c39cafac7b398dad635..b44684a9a696c36f2312ae78ea4018bc261c0d5e 100644
--- a/NonLinearSolver/internalPoints/IntPtData.h
+++ b/NonLinearSolver/internalPoints/IntPtData.h
@@ -55,10 +55,10 @@ class IntPtData{
       return *this;
     };
 		
-		std::vector<ValType>& f(FunctionSpaceBase* space, MElement* ele, IntPt &GP){
+		std::vector<ValType>& f(const FunctionSpaceBase* space, MElement* ele, IntPt &GP){
 			if (!_valEstimated){
 				_valEstimated = true;
-				FunctionSpace<T>* spT = dynamic_cast<FunctionSpace<T>*>(space);
+				const FunctionSpace<T>* spT = dynamic_cast<const FunctionSpace<T>*>(space);
 				if (spT!=NULL){
 					spT->f(ele,GP.pt[0],GP.pt[1],GP.pt[2],val);
 				};
@@ -66,10 +66,10 @@ class IntPtData{
 			return val;
 		};
 		
-		std::vector<GradType>& gradf(FunctionSpaceBase* space, MElement* ele, IntPt &GP){
+		std::vector<GradType>& gradf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP){
 			if (!_gradEstimated){
 				_gradEstimated=true;
-				FunctionSpace<T>* spT = dynamic_cast<FunctionSpace<T>*>(space);
+				const FunctionSpace<T>* spT = dynamic_cast<const FunctionSpace<T>*>(space);
 				if (spT!=NULL){
 					spT->gradf(ele,GP.pt[0],GP.pt[1],GP.pt[2],grad);
 				};
@@ -77,10 +77,10 @@ class IntPtData{
 			return grad;
 		};
 		
-		std::vector<HessType>& hessf(FunctionSpaceBase* space, MElement* ele, IntPt &GP){
+		std::vector<HessType>& hessf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP){
 			if (!_hessEstimated){
 				_hessEstimated = true;
-				FunctionSpace<T>* spT = dynamic_cast<FunctionSpace<T>*>(space);
+				const FunctionSpace<T>* spT = dynamic_cast<const FunctionSpace<T>*>(space);
 				if (spT!=NULL){
 					spT->hessf(ele,GP.pt[0],GP.pt[1],GP.pt[2],hess);
 				};
@@ -88,10 +88,10 @@ class IntPtData{
 			return hess;
 		};
 		
-		std::vector<ThirdDevType>&  thirdDevf(FunctionSpaceBase* space, MElement* ele, IntPt &GP){
+		std::vector<ThirdDevType>&  thirdDevf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP){
 			if (!_thirdEstimated){
 				_thirdEstimated = true;
-				FunctionSpace<T>* spT = dynamic_cast<FunctionSpace<T>*>(space);
+				const FunctionSpace<T>* spT = dynamic_cast<const FunctionSpace<T>*>(space);
 				if (spT!=NULL){
 					spT->thirdDevf(ele,GP.pt[0],GP.pt[1],GP.pt[2],third);
 				};
diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 520074ffe3e768b8bc459780e04cc52d39308c98..cc2444f8953a2a6413287c0a18267a864147d0a7 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -484,7 +484,7 @@ double IPField::computeBulkDeformationEnergy(MElement *ele, const partDomain *do
   double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -531,7 +531,7 @@ double IPField::computeInterfaceDeformationEnergy(MElement *ele, const dgPartDom
   double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -623,7 +623,7 @@ double IPField::computeBulkPlasticEnergy(MElement *ele, const partDomain *dom, c
   double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -669,7 +669,7 @@ double IPField::computeInterfacePlasticEnergy(MElement *ele, const dgPartDomain
   double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -761,7 +761,7 @@ double IPField::computeBulkDamageEnergy(MElement *ele, const partDomain *dom, co
   double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -807,7 +807,7 @@ double IPField::computeInterfaceDamageEnergy(MElement *ele, const dgPartDomain *
   double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -898,7 +898,7 @@ double IPField::computeBulkPathFollowingLocalValue(MElement *ele, const partDoma
   double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -943,7 +943,7 @@ double IPField::computeInterfacePathFollowingLocalValue(MElement *ele, const dgP
 	double ener =0.;
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -1076,7 +1076,7 @@ int IPField::computeFractureEnergy(MElement *ele, const dgPartDomain *dom,double
   AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(ele->getNum());
   // Jacobian values;
   std::vector<double> vdetJ(npts);
-  nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
+  const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(dom->getFunctionSpace());
   if(nlsspace!=NULL)
   {
     std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
@@ -1459,7 +1459,7 @@ double IPField::getIPcompAndVolume(const partDomain *ef,MElement *ele, const IPS
     AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(ele->getNum());
     // Jacobian values;
     std::vector<double> vdetJ(npts);
-    nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<nlsFunctionSpaceUVW<double>*>(ef->getFunctionSpace());
+    const nlsFunctionSpaceUVW<double>* nlsspace = dynamic_cast<const nlsFunctionSpaceUVW<double>*>(ef->getFunctionSpace());
     if(nlsspace!=NULL)
     {
       std::vector<std::vector<TensorialTraits<double>::GradType> > vgradsuvw;
diff --git a/NonLinearSolver/internalPoints/ipFiniteStrain.cpp b/NonLinearSolver/internalPoints/ipFiniteStrain.cpp
index dba79dcc209d37fa368fa63b0ea55b564b29f118..f078b9e0bf1b3803c55969cc63e8fd1cc49bc7d3 100644
--- a/NonLinearSolver/internalPoints/ipFiniteStrain.cpp
+++ b/NonLinearSolver/internalPoints/ipFiniteStrain.cpp
@@ -29,27 +29,27 @@ ipFiniteStrain::~ipFiniteStrain(){
 	}
 }
   // get shape function information from ipvariable
-std::vector<TensorialTraits<double>::ValType>& ipFiniteStrain::f(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::ValType>& ipFiniteStrain::f(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
   if (_GPData == NULL){
 		_GPData = new IntPtData<double>();
 	};
   return _GPData->f(space,ele,GP);
 };
-std::vector<TensorialTraits<double>::GradType>& ipFiniteStrain::gradf(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::GradType>& ipFiniteStrain::gradf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
   if (_GPData == NULL){
 		_GPData = new IntPtData<double>();
 	};
 	return _GPData->gradf(space,ele,GP);
 };
 
-std::vector<TensorialTraits<double>::HessType>& ipFiniteStrain::hessf(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::HessType>& ipFiniteStrain::hessf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
   if (_GPData == NULL){
 		_GPData = new IntPtData<double>();
 	};
 	return _GPData->hessf(space,ele,GP);
 };
 
-std::vector<TensorialTraits<double>::ThirdDevType>& ipFiniteStrain::thirdDevf(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::ThirdDevType>& ipFiniteStrain::thirdDevf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const {
   if (_GPData == NULL){
 		_GPData = new IntPtData<double>();
 	};
@@ -64,14 +64,14 @@ double& ipFiniteStrain::getJacobianDeterminant(MElement* ele, IntPt &GP) const{
 };
 
 
-std::vector<TensorialTraits<double>::ValType>& ipFiniteStrain::f(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::ValType>& ipFiniteStrain::f(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
   if (_interfaceGPData == NULL){
 		_interfaceGPData = new IntPtData<double>();
 	};
 	MElement* e = dynamic_cast<MElement*>(ele);
   return _interfaceGPData->f(space,e,GP);
 };
-std::vector<TensorialTraits<double>::GradType>& ipFiniteStrain::gradf(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::GradType>& ipFiniteStrain::gradf(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
   if (_interfaceGPData == NULL){
 		_interfaceGPData = new IntPtData<double>();
 	};
@@ -79,7 +79,7 @@ std::vector<TensorialTraits<double>::GradType>& ipFiniteStrain::gradf(FunctionSp
 	return _interfaceGPData->gradf(space,e,GP);
 };
 
-std::vector<TensorialTraits<double>::HessType>& ipFiniteStrain::hessf(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::HessType>& ipFiniteStrain::hessf(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
   if (_interfaceGPData == NULL){
 		_interfaceGPData = new IntPtData<double>();
 	};
@@ -87,7 +87,7 @@ std::vector<TensorialTraits<double>::HessType>& ipFiniteStrain::hessf(FunctionSp
 	return _interfaceGPData->hessf(space,e,GP);
 };
 
-std::vector<TensorialTraits<double>::ThirdDevType>& ipFiniteStrain::thirdDevf(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
+std::vector<TensorialTraits<double>::ThirdDevType>& ipFiniteStrain::thirdDevf(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const {
   if (_interfaceGPData == NULL){
 		_interfaceGPData = new IntPtData<double>();
 	};
diff --git a/NonLinearSolver/internalPoints/ipFiniteStrain.h b/NonLinearSolver/internalPoints/ipFiniteStrain.h
index 19ae750e731209ddddaf96e7f42bd1ff40429865..78c943ef5d7b8bbb18aaebb80cf77aedc643a50e 100644
--- a/NonLinearSolver/internalPoints/ipFiniteStrain.h
+++ b/NonLinearSolver/internalPoints/ipFiniteStrain.h
@@ -26,16 +26,16 @@ class ipFiniteStrain : public IPVariableMechanics{
     virtual ~ipFiniteStrain();
 
       // get shape function information from ipvariable
-    virtual std::vector<TensorialTraits<double>::ValType>& f(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
-    virtual std::vector<TensorialTraits<double>::GradType>& gradf(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
-    virtual std::vector<TensorialTraits<double>::HessType>& hessf(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
-    virtual std::vector<TensorialTraits<double>::ThirdDevType>& thirdDevf(FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
+    virtual std::vector<TensorialTraits<double>::ValType>& f(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
+    virtual std::vector<TensorialTraits<double>::GradType>& gradf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
+    virtual std::vector<TensorialTraits<double>::HessType>& hessf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
+    virtual std::vector<TensorialTraits<double>::ThirdDevType>& thirdDevf(const FunctionSpaceBase* space, MElement* ele, IntPt &GP) const;
     virtual double& getJacobianDeterminant(MElement* ele, IntPt &GP) const;
 
-		virtual std::vector<TensorialTraits<double>::ValType>& f(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
-    virtual std::vector<TensorialTraits<double>::GradType>& gradf(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
-    virtual std::vector<TensorialTraits<double>::HessType>& hessf(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
-    virtual std::vector<TensorialTraits<double>::ThirdDevType>& thirdDevf(FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
+		virtual std::vector<TensorialTraits<double>::ValType>& f(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
+    virtual std::vector<TensorialTraits<double>::GradType>& gradf(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
+    virtual std::vector<TensorialTraits<double>::HessType>& hessf(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
+    virtual std::vector<TensorialTraits<double>::ThirdDevType>& thirdDevf(const FunctionSpaceBase* space, MInterfaceElement* ele, IntPt &GP) const;
     virtual double& getJacobianDeterminant(MInterfaceElement* ele, IntPt &GP) const;
 
 
diff --git a/NonLinearSolver/space/ThreeDLagrangeFunctionSpace.h b/NonLinearSolver/space/ThreeDLagrangeFunctionSpace.h
index 1c42c0e769a11e6c33582c3fe19727ee87baf0b1..5e37981e2f000dad26485fbf1a32f4db4a2de3c0 100644
--- a/NonLinearSolver/space/ThreeDLagrangeFunctionSpace.h
+++ b/NonLinearSolver/space/ThreeDLagrangeFunctionSpace.h
@@ -18,16 +18,16 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
 {
  protected :
   std::vector<int> comp;
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys)=0;
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const=0;
   const int _ifield;
  private: // avoid constant reallocation
-   double _ival[256];
-   double _igrad[256][3];
-   double _ihess [256][3][3];
-   double _ithird[256][3][3][3];
-   HessType hesst;
-   GradType gradt;
-   ThirdDevType thirdt;
+   mutable double _ival[256];
+   mutable double _igrad[256][3];
+   mutable double _ihess [256][3][3];
+   mutable double _ithird[256][3][3][3];
+   mutable HessType hesst;
+   mutable GradType gradt;
+   mutable ThirdDevType thirdt;
  public:
 
   ThreeDLagrangeFunctionSpace(int id, int ncomp,const bool hesscompute,
@@ -58,7 +58,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
   }
   virtual ~ThreeDLagrangeFunctionSpace(){}
 
-  virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+  virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) const
   {
     ele->getShapeFunctions(u,v,w,_ival);
     int nbFF = ele->getNumShapeFunctions();
@@ -69,7 +69,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
         vals[valssize + i+l*nbFF] = _ival[i];
   }
 
-  virtual void fuvw(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+  virtual void fuvw(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) const
   {
     ele->getShapeFunctions(u,v,w,_ival);
     int nbFF = ele->getNumShapeFunctions();
@@ -80,7 +80,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
         vals[valssize + i+l*nbFF] = _ival[i];
   }
 
-  virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+  virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) const
   {
     ele->getGradShapeFunctions(u, v, w, _igrad);
     int nbFF = ele->getNumShapeFunctions();
@@ -99,7 +99,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
         grads[gradssize + i+l*nbFF]=gradt;
     }
   }
-  virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
+  virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) const
   {
     ele->getGradShapeFunctions(u,v,w,_igrad);
     int nbFF = ele->getNumShapeFunctions();
@@ -113,7 +113,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
         grads[gradssize + i + l*nbFF] = gradt;
     }
   }
-  virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+  virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess) const
   {
     ele->getHessShapeFunctions(u,v,w,_ihess);
     int nbFF = ele->getNumShapeFunctions();
@@ -128,7 +128,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
       }
     }
   }
-  virtual void hessf(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+  virtual void hessf(MElement *ele, double u, double v, double w,std::vector<HessType> &hess) const
   {
     ele->getGradShapeFunctions(u,v,w,_igrad);
     ele->getHessShapeFunctions(u,v,w,_ihess);
@@ -177,7 +177,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
     };
   };
 
-  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,std::vector<ThirdDevType> &third){
+  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,std::vector<ThirdDevType> &third) const{
     ele->getThirdDerivativeShapeFunctions(u,v,w,_ithird);
     int nbFF = ele->getNumShapeFunctions();
     int thirdsize = third.size();
@@ -192,7 +192,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
       }
     }
   }; //need to high order fem
-  virtual void thirdDevf(MElement *ele, double u, double v, double w,std::vector<ThirdDevType> &third){
+  virtual void thirdDevf(MElement *ele, double u, double v, double w,std::vector<ThirdDevType> &third) const{
     ele->getGradShapeFunctions(u,v,w,_igrad);
     ele->getHessShapeFunctions(u,v,w,_ihess);
     ele->getThirdDerivativeShapeFunctions(u,v,w,_ithird);
@@ -270,10 +270,10 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
     };
   }; //need to high order fem
 
-  virtual int getNumKeys(MElement *ele) {if (ele->getParent()) ele = ele->getParent();return _ncomp*ele->getNumVertices();}
+  virtual int getNumKeys(MElement *ele)  const {if (ele->getParent()) ele = ele->getParent();return _ncomp*ele->getNumVertices();}
   virtual int getId(void) const {return _ifield;}
 
-  virtual void getKeys(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const{
   /*
     // Keys depends if MElement *ele is a bulk element or an interface element
     if(ele->getDim() == 2){
@@ -303,7 +303,7 @@ class ThreeDLagrangeFunctionSpace : public nlsFunctionSpaceXYZ<double>
 	virtual FunctionSpaceBase* clone(const int id) const = 0;
 
  protected :
-  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys){
+  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys) const{
     this->getKeysOnElement(iele->getElem(0),keys);
     if(!(iele->getElem(0) == iele->getElem(1)))
       this->getKeysOnElement(iele->getElem(1),keys);
@@ -316,15 +316,15 @@ class IsoparametricLagrangeFunctionSpace : public nlsFunctionSpaceUVW<double>
 {
  protected :
   std::vector<int> comp;
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) {Msg::Error("You forget to define getKeysOnElement for your functional space");}
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const {Msg::Error("You forget to define getKeysOnElement for your functional space");}
   const int _ifield;
  public:
  private: // avoid constant reallocation
-   double _ival[256];
-   double _igrad[256][3];
-   double _ihess [256][3][3];
-   HessType hesst;
-   GradType gradt;
+   mutable double _ival[256];
+   mutable double _igrad[256][3];
+   mutable double _ihess [256][3][3];
+   mutable HessType hesst;
+   mutable GradType gradt;
  public:
 
   IsoparametricLagrangeFunctionSpace(int id) : nlsFunctionSpaceUVW<double>(3,true,false), _ifield(id)
@@ -354,7 +354,7 @@ class IsoparametricLagrangeFunctionSpace : public nlsFunctionSpaceUVW<double>
 		}
 	};
 
-  virtual void fuvw(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
+  virtual void fuvw(MElement *ele, double u, double v, double w,std::vector<ValType> &vals) const
   {
     ele->getShapeFunctions(u,v,w,_ival);
     int nbFF = ele->getNumShapeFunctions();
@@ -365,12 +365,12 @@ class IsoparametricLagrangeFunctionSpace : public nlsFunctionSpaceUVW<double>
         vals[valssize + i+l*nbFF] = _ival[i];
   }
   // ==fuvw
-  virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+  virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) const
   {
     this->fuvw(ele,u,v,w,vals);
   }
 
-  virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
+  virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) const
   {
     ele->getGradShapeFunctions(u,v,w,_igrad);
     int nbFF = ele->getNumShapeFunctions();
@@ -385,12 +385,12 @@ class IsoparametricLagrangeFunctionSpace : public nlsFunctionSpaceUVW<double>
     }
   }
   // == gradfuvw
-  virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+  virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) const
   {
     this->gradfuvw(ele,u,v,w,grads);
   }
 
-  virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+  virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess) const
   {
     ele->getHessShapeFunctions(u,v,w,_ihess);
     int nbFF = ele->getNumShapeFunctions();
@@ -406,20 +406,20 @@ class IsoparametricLagrangeFunctionSpace : public nlsFunctionSpaceUVW<double>
     }
   }
   // == hessfuvw
-  virtual void hessf(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+  virtual void hessf(MElement *ele, double u, double v, double w,std::vector<HessType> &hess) const
   {
     this->hessfuvw(ele,u,v,w,hess);
   }
-  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,ContainerThirdDevType &third){}; // Empty as "false"
+  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,ContainerThirdDevType &third) const{}; // Empty as "false"
   // == thirdDevuvw
-  virtual void thirdDevf(MElement *ele, double u, double v, double w,ContainerThirdDevType &third)
+  virtual void thirdDevf(MElement *ele, double u, double v, double w,ContainerThirdDevType &third) const
   {
     this->thirdDevfuvw(ele,u,v,w,third);
   }
-  virtual int getNumKeys(MElement *ele) {if (ele->getParent()) ele = ele->getParent();return _ncomp*ele->getNumVertices();}
+  virtual int getNumKeys(MElement *ele) const {if (ele->getParent()) ele = ele->getParent();return _ncomp*ele->getNumVertices();}
   virtual int getId(void) const {return _ifield;}
 
-  virtual void getKeys(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const{
     // Keys depends if MElement *ele is a bulk element or an interface element
     if(ele->getDim() == 1){
       MInterfaceElement *iele = dynamic_cast<MInterfaceElement*>(ele);
@@ -441,7 +441,7 @@ class IsoparametricLagrangeFunctionSpace : public nlsFunctionSpaceUVW<double>
   }
 
  protected :
-  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys){
+  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys) const{
     this->getKeysOnElement(iele->getElem(0),keys);
     if(!(iele->getElem(0) == iele->getElem(1)))
       this->getKeysOnElement(iele->getElem(1),keys);
diff --git a/NonLinearSolver/space/interFunctionSpace.h b/NonLinearSolver/space/interFunctionSpace.h
index 97849eb0c6f708dfdf73394ccc15abddbcee4563..e0cda6d2cf06960cdd462234de94f17b8fc42025 100644
--- a/NonLinearSolver/space/interFunctionSpace.h
+++ b/NonLinearSolver/space/interFunctionSpace.h
@@ -13,11 +13,13 @@
 #include "functionSpace.h"
 class interFunctionSpace
 {
- public:
-  virtual FunctionSpaceBase* getMinusSpace() const=0;
-  virtual FunctionSpaceBase* getPlusSpace() const=0;
-  virtual void getNumKeys(MInterfaceElement *ele, int &numMinus, int &numPlus) const= 0; // if one needs the number of dofs
-  virtual void getKeys(MInterfaceElement *ele, std::vector<Dof> &Rminus,std::vector<Dof> &Rplus) const=0;
+  public:
+    virtual FunctionSpaceBase* getMinusSpace() =0;
+    virtual FunctionSpaceBase* getPlusSpace() =0;
+    virtual const FunctionSpaceBase* getMinusSpace() const=0;
+    virtual const FunctionSpaceBase* getPlusSpace() const=0;
+    virtual void getNumKeys(MInterfaceElement *ele, int &numMinus, int &numPlus) const= 0; // if one needs the number of dofs
+    virtual void getKeys(MInterfaceElement *ele, std::vector<Dof> &Rminus,std::vector<Dof> &Rplus) const=0;
 };
 
 #endif //_INTERFUNCTIONSPACE_H_
diff --git a/NonLinearSolver/space/nlsFunctionSpace.h b/NonLinearSolver/space/nlsFunctionSpace.h
index 498b9de6a6cdf144fbb8850746bcd267eec3ccd2..072b8eca214f284599f17a9c5129a02e3ed17ff4 100644
--- a/NonLinearSolver/space/nlsFunctionSpace.h
+++ b/NonLinearSolver/space/nlsFunctionSpace.h
@@ -46,7 +46,7 @@ template<class T> struct GaussPointSpaceValues
   ContainerHessType _vhess;
   ContainerThirdDevType _vthird;
   double _detJ; // value of detJacobian only for XYZ space as in UVW it is the same container for all elements !!
-  GaussPointSpaceValues(nlsFunctionSpaceUVW<T> *sp,MElement *ele,const int ncomp,const double u, const double v, const double w,
+  GaussPointSpaceValues(const nlsFunctionSpaceUVW<T> *sp,MElement *ele,const int ncomp,const double u, const double v, const double w,
                         const bool hessian, const bool thirdDev) : _detJ(0.)
   {
     // compute the value
@@ -97,10 +97,10 @@ template<class T> class nlsFunctionSpaceUVW : public FunctionSpace<T>
   const bool _hessianComputation;
   const bool _thirdDevComputation;
   // map to store the data
-  std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> > _mapGaussValues;
+  mutable std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> > _mapGaussValues; // cache in order to assess everywhere
 
   // function to allocate the map
-  typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator allocatePoints(MElement *ele,int npts, IntPt *GP)
+  typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator allocatePoints(MElement *ele,int npts, IntPt *GP) const
   {
     std::vector<GaussPointSpaceValues<T>*> allval;
     for(int i=0;i<npts;i++)
@@ -112,12 +112,12 @@ template<class T> class nlsFunctionSpaceUVW : public FunctionSpace<T>
     return pib.first;
   }
  public:
-  virtual int getNumKeys(MElement *ele) = 0;
-  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) = 0;
+  virtual int getNumKeys(MElement *ele) const = 0;
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const = 0;
   virtual void getComp(std::vector<int> &comp) const = 0;
 	virtual FunctionSpaceBase* clone(const int id) const  = 0;
 
-  virtual void getKeysOnVertex(MElement* ele, MVertex* v, const std::vector<int>& vercomp, std::vector<Dof>& keys){
+  virtual void getKeysOnVertex(MElement* ele, MVertex* v, const std::vector<int>& vercomp, std::vector<Dof>& keys) const{
     // get keys on element
     std::vector<Dof> elekeys;
     this->getKeys(ele,elekeys);
@@ -174,19 +174,19 @@ template<class T> class nlsFunctionSpaceUVW : public FunctionSpace<T>
     _mapGaussValues.clear();
   }
   // traditionel functions
-  virtual void fuvw(MElement *ele, double u, double v, double w, ContainerValType &vals)=0;
-  virtual void f(MElement *ele, double u, double v, double w, ContainerValType &vals) = 0;
-  virtual void gradfuvw(MElement *ele, double u, double v, double w, ContainerGradType &grads)=0;
-  virtual void gradf(MElement *ele, double u, double v, double w, ContainerGradType &grads) = 0;
-  virtual void hessfuvw(MElement *ele, double u, double v, double w, ContainerHessType &hess)=0;
-  virtual void hessf(MElement *ele, double u, double v, double w, ContainerHessType &hess) =0; // define it empty if _hessianComputation == false
-  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,ContainerThirdDevType &third)=0;
-  virtual void thirdDevf(MElement *ele, double u, double v, double w,ContainerThirdDevType &third)=0; // define it empty if _thirdDevComputation == false
+  virtual void fuvw(MElement *ele, double u, double v, double w, ContainerValType &vals) const=0;
+  virtual void f(MElement *ele, double u, double v, double w, ContainerValType &vals) const= 0;
+  virtual void gradfuvw(MElement *ele, double u, double v, double w, ContainerGradType &grads) const=0;
+  virtual void gradf(MElement *ele, double u, double v, double w, ContainerGradType &grads) const = 0;
+  virtual void hessfuvw(MElement *ele, double u, double v, double w, ContainerHessType &hess) const=0;
+  virtual void hessf(MElement *ele, double u, double v, double w, ContainerHessType &hess) const =0; // define it empty if _hessianComputation == false
+  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,ContainerThirdDevType &third) const=0;
+  virtual void thirdDevf(MElement *ele, double u, double v, double w,ContainerThirdDevType &third) const=0; // define it empty if _thirdDevComputation == false
 
 
   // Functions needed to use a map which store the values for the different Gauss Points
   // These ones can be const if fuvw,gradfuvw, hessfuvw and thirdfuvw are const functions
-  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< ContainerValType > &vvals)
+  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< ContainerValType > &vvals) const
   {
     typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator itG = _mapGaussValues.find(GP);
     // fill the value if not allocated
@@ -198,7 +198,7 @@ template<class T> class nlsFunctionSpaceUVW : public FunctionSpace<T>
     for(int i=0;i<npts;i++)
       vvals.push_back(allc[i]->_vvals); // space in uvw --> vvals == _vvals;
   }
-  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< GaussPointSpaceValues<T>*> &vall)
+  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< GaussPointSpaceValues<T>*> &vall) const
   {
     typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator itG = _mapGaussValues.find(GP);
     // fill the value if not allocated
@@ -209,7 +209,7 @@ template<class T> class nlsFunctionSpaceUVW : public FunctionSpace<T>
     std::vector<GaussPointSpaceValues<T>*>& allc = itG->second;
     vall.assign(allc.begin(),allc.end()); // space in uvw --> vall == allc;
   }
-  virtual void getgradf(MElement *ele, int npts, IntPt *GP, std::vector< ContainerGradType > &vgrads)
+  virtual void getgradf(MElement *ele, int npts, IntPt *GP, std::vector< ContainerGradType > &vgrads) const
   {
     typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator itG = _mapGaussValues.find(GP);
     // fill the value if not allocated
@@ -222,11 +222,11 @@ template<class T> class nlsFunctionSpaceUVW : public FunctionSpace<T>
       vgrads.push_back(allc[i]->_vgrads); // space in uvw --> vgrads == _vgrads;
   }
   // Needed to compute the Jacobian quickly !
-  virtual void getgradfuvw(MElement *ele, int npts, IntPt *GP, std::vector< ContainerGradType > &vgrads)
+  virtual void getgradfuvw(MElement *ele, int npts, IntPt *GP, std::vector< ContainerGradType > &vgrads) const
   {
     this->getgradf(ele,npts,GP,vgrads);
   }
-  virtual void gethessf(MElement *ele, int npts, IntPt *GP, std::vector< ContainerHessType > &vhesss)
+  virtual void gethessf(MElement *ele, int npts, IntPt *GP, std::vector< ContainerHessType > &vhesss) const
   {
     if(_hessianComputation)
     {
@@ -244,7 +244,7 @@ template<class T> class nlsFunctionSpaceUVW : public FunctionSpace<T>
       Msg::Fatal("You compute the hessian but give false for _hessianComputation to your function space");
     }
   }
-  virtual void getthirdDevf(MElement *ele, int npts, IntPt *GP,std::vector<ContainerThirdDevType> &third)
+  virtual void getthirdDevf(MElement *ele, int npts, IntPt *GP,std::vector<ContainerThirdDevType> &third) const
   {
     if(_thirdDevComputation)
     {
@@ -280,9 +280,9 @@ template<class T> class nlsFunctionSpaceXYZ : public nlsFunctionSpaceUVW<T>
   mutable GradType gradt;
   mutable HessType hesst;
   mutable ThirdDevType thirdt;
-  std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> > _mapgt; // to avoid constant allocation FIX double map search (one uvw one xyz) HOW ??
+  mutable std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> > _mapgt; // to avoid constant allocation FIX double map search (one uvw one xyz) HOW ??
   // allocation of the map
-  typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator getStorageSpace(MElement *ele,const int npts,IntPt* GP)
+  typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator getStorageSpace(MElement *ele,const int npts,IntPt* GP) const
   {
     typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator itG = _mapgt.find(GP);
     if(itG == _mapgt.end())
@@ -299,8 +299,8 @@ template<class T> class nlsFunctionSpaceXYZ : public nlsFunctionSpaceUVW<T>
     return itG;
   }
  public:
-  virtual int getNumKeys(MElement *ele) = 0;
-  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) = 0;
+  virtual int getNumKeys(MElement *ele) const= 0;
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const = 0;
   virtual void getComp(std::vector<int> &comp) const = 0;
   nlsFunctionSpaceXYZ(const int nc,const bool hessian,const bool thirdDev) : nlsFunctionSpaceUVW<T>(nc,hessian,thirdDev),
                                                                              gradt(0.), hesst(0.), thirdt(0.){}
@@ -316,17 +316,17 @@ template<class T> class nlsFunctionSpaceXYZ : public nlsFunctionSpaceUVW<T>
     _mapgt.clear();
   }
   // traditionel functions
-  virtual void fuvw(MElement *ele, double u, double v, double w, ContainerValType &vals)=0;
-  virtual void f(MElement *ele, double u, double v, double w, ContainerValType &vals) = 0;
-  virtual void gradfuvw(MElement *ele, double u, double v, double w, ContainerGradType &grads)=0;
-  virtual void gradf(MElement *ele, double u, double v, double w, ContainerGradType &grads) = 0;
-  virtual void hessfuvw(MElement *ele, double u, double v, double w, ContainerHessType &hess)=0;
-  virtual void hessf(MElement *ele, double u, double v, double w, ContainerHessType &hess) =0; // define it empty if _hessianComputation == false
-  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,ContainerThirdDevType &third)=0;
-  virtual void thirdDevf(MElement *ele, double u, double v, double w,ContainerThirdDevType &third)=0; // define it empty if _thirdDevComputation == false
+  virtual void fuvw(MElement *ele, double u, double v, double w, ContainerValType &vals) const=0;
+  virtual void f(MElement *ele, double u, double v, double w, ContainerValType &vals) const= 0;
+  virtual void gradfuvw(MElement *ele, double u, double v, double w, ContainerGradType &grads)const=0;
+  virtual void gradf(MElement *ele, double u, double v, double w, ContainerGradType &grads) const= 0;
+  virtual void hessfuvw(MElement *ele, double u, double v, double w, ContainerHessType &hess) const=0;
+  virtual void hessf(MElement *ele, double u, double v, double w, ContainerHessType &hess) const =0; // define it empty if _hessianComputation == false
+  virtual void thirdDevfuvw(MElement *ele, double u, double v, double w,ContainerThirdDevType &third) const=0;
+  virtual void thirdDevf(MElement *ele, double u, double v, double w,ContainerThirdDevType &third) const=0; // define it empty if _thirdDevComputation == false
 
   // Redifine the functions to return the values in XYZ
-  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< ContainerValType > &vvals)
+  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< ContainerValType > &vvals) const
   {
     typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator itG = this->_mapGaussValues.find(GP);
     // fill the value if not allocated
@@ -340,7 +340,7 @@ template<class T> class nlsFunctionSpaceXYZ : public nlsFunctionSpaceUVW<T>
       vvals.push_back(allc[i]->_vvals); // value in uvw == value in xyz
   }
   //  no resize allowed after initialization --> grad.resize() is not allowed idem val,hess and third !!
-  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< GaussPointSpaceValues<T>*> &vall)
+  virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< GaussPointSpaceValues<T>*> &vall) const
   {
     typename std::map<IntPt*,std::vector<GaussPointSpaceValues<T>*> >::iterator itG = this->_mapGaussValues.find(GP);
     // fill the value if not allocated
@@ -483,7 +483,7 @@ template<class T> class nlsFunctionSpaceXYZ : public nlsFunctionSpaceUVW<T>
     vall.assign(allc.begin(),allc.end()); // vall have pointer on storage space
   }
   // Needed to compute the Jacobian quickly without computing third
-  virtual void getgradfuvw(MElement *ele, int npts, IntPt *GP, std::vector< ContainerGradType > &vgrads)
+  virtual void getgradfuvw(MElement *ele, int npts, IntPt *GP, std::vector< ContainerGradType > &vgrads) const
   {
     nlsFunctionSpaceUVW<T>::getgradf(ele,npts,GP,vgrads);
   }
@@ -521,7 +521,7 @@ template<class T> class nlsFunctionSpaceXYZ : public nlsFunctionSpaceUVW<T>
       }
     }
   }
-  virtual void gethessf(MElement *ele, int npts, IntPt *GP, std::vector< ContainerHessType > &vhesss)
+  virtual void gethessf(MElement *ele, int npts, IntPt *GP, std::vector< ContainerHessType > &vhesss) const
   {
     if(this->_hessianComputation)
     {
@@ -590,7 +590,7 @@ template<class T> class nlsFunctionSpaceXYZ : public nlsFunctionSpaceUVW<T>
       Msg::Fatal("You compute the hessian but give false for _hessianComputation to your function space");
     }
   }
-  virtual void getthirdDevf(MElement *ele, int npts, IntPt *GP,std::vector<ContainerThirdDevType> &vthird)
+  virtual void getthirdDevf(MElement *ele, int npts, IntPt *GP,std::vector<ContainerThirdDevType> &vthird) const
   {
     if(this->_thirdDevComputation)
     {
diff --git a/dG3D/src/dG3DContactFunctionSpace.h b/dG3D/src/dG3DContactFunctionSpace.h
index fcfb9856d343ff549506db70707e5ded36fec2fc..331c7695ffc2f156d39e4573f570fe484e4f3fee 100644
--- a/dG3D/src/dG3DContactFunctionSpace.h
+++ b/dG3D/src/dG3DContactFunctionSpace.h
@@ -24,7 +24,7 @@ class dG3DRigidContactSpace : public rigidContactSpaceBase{
   dG3DRigidContactSpace(const int id, FunctionSpace<double> *sp, MVertex *ver, const int comp1,
                    const int comp2 = -1, const int comp3 =-1) : rigidContactSpaceBase(id,sp,ver,comp1,comp2,comp3){}
   virtual ~dG3DRigidContactSpace(){}
-  virtual void getKeys(MElement *ele,std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele,std::vector<Dof> &keys) const{
     if(ele->getDim() > 0){ // if dim 0 return the key of gravity center only !!
       _space->getKeys(ele,keys);
     }
@@ -32,18 +32,18 @@ class dG3DRigidContactSpace : public rigidContactSpaceBase{
       keys.push_back(Dof(_gc->getNum(),dG3DDof3IntType::createTypeWithThreeInts(_comp[i],_id)));
     }
   }
-  virtual int getNumKeys(MElement *ele){
+  virtual int getNumKeys(MElement *ele) const{
     int nkeysele = _space->getNumKeys(ele);
     return nkeysele + _comp.size();
   }
-  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys){
+  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys) const{
     for(int i=0;i<_comp.size(); i++)
       keys.push_back(Dof(_gc->getNum(),dG3DDof3IntType::createTypeWithThreeInts(_comp[i],_id)));
   }
-  virtual int getNumKeysOfGravityCenter(){
+  virtual int getNumKeysOfGravityCenter() const{
     return _comp.size();
   }
-  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys) const{
     // generate keys of element and select the good ones after LAGRANGE OK ?? CHANGE THIS HOW TODO
     std::vector<Dof> tkeys;
     this->getKeys(ele,tkeys);
@@ -54,7 +54,7 @@ class dG3DRigidContactSpace : public rigidContactSpaceBase{
       keys.push_back(tkeys[ind+i*nbver]);
     this->getKeysOfGravityCenter(keys);
   }
-  virtual int getNumKeys(MElement *ele, int ind){
+  virtual int getNumKeys(MElement *ele, int ind) const{
     return 2*_comp.size();
   }
 };
diff --git a/dG3D/src/dG3DDomain.h b/dG3D/src/dG3DDomain.h
index 0d2e1bfe9d09670bfc7343afa0ff776065efe85c..4fde36a4e7e0a793958ae20d1560e0e4e88c6e92 100644
--- a/dG3D/src/dG3DDomain.h
+++ b/dG3D/src/dG3DDomain.h
@@ -133,10 +133,14 @@ class dG3DDomain : public dgPartDomain{
                                                     const simpleFunctionTime<double>* f,const unknownField *uf,
                                                     const IPField * ip, nonLinearBoundaryCondition::location onWhat) const;
   virtual double scaleTimeStep() const {return _sts;}
-  virtual FunctionSpaceBase* getFunctionSpace() const{return _space;}
-  virtual FunctionSpaceBase* getFunctionSpaceMinus() const{return _space;}
-  virtual FunctionSpaceBase* getFunctionSpacePlus() const{return _space;}
-  virtual FunctionSpaceBase* getInterfaceFunctionSpace() const {return _extraSpace;};
+  virtual FunctionSpaceBase* getFunctionSpace(){return _space;}
+  virtual const FunctionSpaceBase* getFunctionSpace() const {return _space;}
+  virtual FunctionSpaceBase* getFunctionSpaceMinus() {return _space;}
+  virtual const FunctionSpaceBase* getFunctionSpaceMinus() const{return _space;}
+  virtual FunctionSpaceBase* getFunctionSpacePlus() {return _space;}
+  virtual const FunctionSpaceBase* getFunctionSpacePlus() const{return _space;}
+  virtual FunctionSpaceBase* getInterfaceFunctionSpace() {return _extraSpace;};
+  virtual const FunctionSpaceBase* getInterfaceFunctionSpace() const {return _extraSpace;};
   virtual const partDomain* getMinusDomain() const{return this;}
   virtual const partDomain* getPlusDomain() const{return this;}
   virtual partDomain* getMinusDomain() {return this;}
@@ -275,8 +279,10 @@ class interDomainBetween3D : public interDomainBase, public dG3DDomain{
   virtual int getPlusLawNum() const{return _domPlus->getLawNum();}
   virtual const partDomain* getMinusDomain() const{return _domMinus;}
   virtual const partDomain* getPlusDomain() const{return _domPlus;}
-	virtual FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
-  virtual FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
+	virtual FunctionSpaceBase* getFunctionSpaceMinus(){return _spaceMinus;}
+  virtual const FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
+  virtual FunctionSpaceBase* getFunctionSpacePlus(){return _spacePlus;}
+  virtual const FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
   virtual partDomain* getMinusDomain(){return _domMinus;}
   virtual partDomain* getPlusDomain() {return _domPlus;}
   void initializeTerms(unknownField *uf,IPField *ip);
@@ -319,8 +325,10 @@ class nonLocalInterDomainBetween3D : public interDomainBase, public nonLocalDama
   virtual const partDomain* getPlusDomain() const{return _domPlus;}
   virtual partDomain* getMinusDomain(){return _domMinus;}
   virtual partDomain* getPlusDomain() {return _domPlus;}
-	virtual FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
-  virtual FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
+	virtual FunctionSpaceBase* getFunctionSpaceMinus(){return _spaceMinus;}
+  virtual const FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
+  virtual FunctionSpaceBase* getFunctionSpacePlus(){return _spacePlus;}
+  virtual const FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
   void initializeTerms(unknownField *uf,IPField *ip);
 
   virtual void createTerms(unknownField *uf,IPField*ip);
@@ -424,8 +432,10 @@ class ExtraDofDiffusionInterDomainBetween3D : public interDomainBase, public Ext
   virtual const partDomain* getPlusDomain() const{return _domPlus;}
   virtual partDomain* getMinusDomain(){return _domMinus;}
   virtual partDomain* getPlusDomain() {return _domPlus;}
-	virtual FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
-  virtual FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
+	virtual FunctionSpaceBase* getFunctionSpaceMinus() {return _spaceMinus;}
+  virtual const FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
+  virtual FunctionSpaceBase* getFunctionSpacePlus() {return _spacePlus;}
+  virtual const FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
   void initializeTerms(unknownField *uf,IPField *ip);
 
   virtual void createTerms(unknownField *uf,IPField*ip)=0;
@@ -597,8 +607,10 @@ class ElecTherMechInterDomainBetween3D : public interDomainBase, public ElecTher
   virtual const partDomain* getPlusDomain() const{return _domPlus;}
   virtual partDomain* getMinusDomain(){return _domMinus;}
   virtual partDomain* getPlusDomain() {return _domPlus;}
-	virtual FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
-  virtual FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
+	virtual FunctionSpaceBase* getFunctionSpaceMinus() {return _spaceMinus;}
+  virtual const FunctionSpaceBase* getFunctionSpaceMinus()  const {return _spaceMinus;}
+  virtual FunctionSpaceBase* getFunctionSpacePlus(){return _spacePlus;}
+  virtual const FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
   void initializeTerms(unknownField *uf,IPField *ip);
 
   virtual void createTerms(unknownField *uf,IPField*ip);
diff --git a/dG3D/src/dG3DFunctionSpace.h b/dG3D/src/dG3DFunctionSpace.h
index 98b8019ee747b874f78727998ae53cc061577bcc..a44b5d3ce43cafb14c0b2b74f921dbfa71d3cded 100644
--- a/dG3D/src/dG3DFunctionSpace.h
+++ b/dG3D/src/dG3DFunctionSpace.h
@@ -33,7 +33,7 @@ class g3DLagrangeFunctionSpace : public ThreeDLagrangeFunctionSpace{
   virtual ~g3DLagrangeFunctionSpace(){};
 
  protected :
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys)
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const
   {
     int nk=ele->getNumVertices(); // return the number of vertices
     // negative type in mpi if the Dof are not located on this partition
@@ -104,7 +104,7 @@ class g3DLagrangeFunctionSpace : public ThreeDLagrangeFunctionSpace{
 
 class dG3DLagrangeFunctionSpace : public ThreeDLagrangeFunctionSpace{
  protected:
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const{
     int nk=ele->getNumVertices(); // return the number of vertices
     #if defined(HAVE_MPI) // small duplication to avoid multiple if in a loop
     if( (ele->getPartition() != 0)and (ele->getPartition() != Msg::GetCommRank() +1))
@@ -169,7 +169,7 @@ class g3DhoDGFunctionSpace : public ThreeDLagrangeFunctionSpace{
   g3DhoDGFunctionSpace(int id, int ncomp, int comp1, int comp2, int comp3) : ThreeDLagrangeFunctionSpace(id,ncomp,comp1,comp2,comp3,true,true){};
   virtual ~g3DhoDGFunctionSpace(){};
  protected :
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys)
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const
   {
     int nk=ele->getNumVertices(); // return the number of vertices
     // negative type in mpi if the Dof are not located on this partition
@@ -211,7 +211,7 @@ class dG3DhoDGFunctionSpace : public ThreeDLagrangeFunctionSpace{
   virtual ~dG3DhoDGFunctionSpace(){};
 
  protected: // avoid duplication of the following functions with dG3DLagrangeFunctionSpace HOW ??
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const{
     int nk=ele->getNumVertices(); // return the number of vertices
     #if defined(HAVE_MPI) // small duplication to avoid multiple if in a loop
     if( (ele->getPartition() != 0)and (ele->getPartition() != Msg::GetCommRank() +1))
@@ -263,7 +263,7 @@ class g3DDirichletBoundaryConditionLagrangeFunctionSpace : public g3DLagrangeFun
 		else return new g3DDirichletBoundaryConditionLagrangeFunctionSpace(id,_ncomp);
  	}
 
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys)
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const
   {
     int nk=ele->getNumVertices();
     // negative type in mpi if the Dof are not located on this partition
@@ -349,7 +349,7 @@ class g3DNeumannBoundaryConditionLagrangeFunctionSpace : public g3DLagrangeFunct
 		else return new g3DNeumannBoundaryConditionLagrangeFunctionSpace(id,_ncomp);
 	}
 
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys)
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const
   {
      // For Neumann BC we sure that the BC is applied on element that are on this partition.
      // Indeed Neumann BC are never applied on ghost element
@@ -381,7 +381,7 @@ class g3DDirichletBoundaryConditionLagrangeFunctionSpaceGhost : public g3DDirich
 	} 
 
   // be sure that the type is negative !!
-  virtual void getKeys(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const{
     int nk=ele->getNumVertices();
     if(ele->getPartition()!=0){ // Sure than it will match
       for (int j=0;j<_ncomp;++j)
@@ -1509,15 +1509,15 @@ class dG3DBoundaryConditionLagrangeFunctionSpace : public dG3DLagrangeFunctionSp
 		Msg::Fatal("dG3DBoundaryConditionLagrangeFunctionSpace::clone needs to be defined");
 	}
 
-  virtual void getKeys(MInterfaceElement *ielem, std::vector<Dof> &keys){
+  virtual void getKeys(MInterfaceElement *ielem, std::vector<Dof> &keys) const{
     Msg::Error("Impossible to get keys on interface element for a Dirichlet Boundary Conditions");
   }
-  virtual void getKeys(MElement *e, std::vector<Dof> &keys)
+  virtual void getKeys(MElement *e, std::vector<Dof> &keys) const
   {
-    std::map<MElement*,std::pair<MElement*,MElement*> >::iterator itele=mapInter.find(e);
+    std::map<MElement*,std::pair<MElement*,MElement*> >::const_iterator itele=mapInter.find(e);
     MElement *ele1 = itele->second.first;
     MElement *ele2 = itele->second.second;
-    std::map<MElement*,std::pair<std::vector<int>,std::vector<int> > >::iterator itvec=mapLocalVertex.find(e);
+    std::map<MElement*,std::pair<std::vector<int>,std::vector<int> > >::const_iterator itvec=mapLocalVertex.find(e);
     int nbcomp=comp.size();
     #if defined(HAVE_MPI) // small duplication to avoid multiple if in a loop
     if( (ele1->getPartition() != 0)and (ele1->getPartition() != Msg::GetCommRank() +1))
@@ -1566,7 +1566,7 @@ class dG3DLagrangeBetween2DomainsFunctionSpace : public ThreeDLagrangeFunctionSp
  protected:
   FunctionSpaceBase *spaceMinus;
   FunctionSpaceBase *spacePlus;
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const{
     Msg::Error("Impossible to get key for an element on a interface Domain");
   }
  public:
@@ -1578,18 +1578,23 @@ class dG3DLagrangeBetween2DomainsFunctionSpace : public ThreeDLagrangeFunctionSp
 		return new dG3DLagrangeBetween2DomainsFunctionSpace(_ncomp,spaceMinus,spacePlus);
 	}
 
-  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys){
+  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys) const{
     spaceMinus->getKeys(iele->getElem(0),keys);
     spacePlus->getKeys(iele->getElem(1),keys);
   }
-  virtual void getKeys(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const{
     // As all element are interface element on an interface domain
     MInterfaceElement *iele = dynamic_cast<MInterfaceElement*>(ele);
     this->getKeys(iele,keys);
   }
   // special function of interFunctionSpace
-  virtual FunctionSpaceBase* getMinusSpace()const {return spaceMinus;}
-  virtual FunctionSpaceBase* getPlusSpace() const {return spacePlus;}
+  virtual FunctionSpaceBase* getMinusSpace() {return spaceMinus;}
+  virtual FunctionSpaceBase* getPlusSpace()  {return spacePlus;}
+  
+  virtual const FunctionSpaceBase* getMinusSpace()const {return spaceMinus;}
+  virtual const FunctionSpaceBase* getPlusSpace() const {return spacePlus;}
+  
+  
   virtual void getNumKeys(MInterfaceElement *ele, int &numMinus, int &numPlus) const
   {
     numMinus = spaceMinus->getNumKeys(ele->getElem(0));
diff --git a/dG3D/src/dG3DHomogenizedTangentTerms.cpp b/dG3D/src/dG3DHomogenizedTangentTerms.cpp
index f20bcadef170e195c260d448548d7b680357c279..4e2aee0983727da4acf18a9c607e9abba339f587 100644
--- a/dG3D/src/dG3DHomogenizedTangentTerms.cpp
+++ b/dG3D/src/dG3DHomogenizedTangentTerms.cpp
@@ -13,7 +13,7 @@
 #include "nonLinearMechSolver.h"
 
 void dG3DHomogenizedTangent::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) const{
-	FunctionSpaceBase* sp = _dom->getFunctionSpace();
+	const FunctionSpaceBase* sp = _dom->getFunctionSpace();
 	
 	int nbDof = sp->getNumKeys(ele);
 	int nbFF = ele->getNumShapeFunctions();
@@ -98,7 +98,7 @@ void dG3DNonLocalHomogenizedTangent::get(MElement *ele,int npts,IntPt *GP,fullMa
 	// get mechanical parts
 	dG3DHomogenizedTangent::get(ele,npts,GP,m);
 	
-	FunctionSpaceBase* space = _dom->getFunctionSpace();
+	const FunctionSpaceBase* space = _dom->getFunctionSpace();
 	
 	int nbDof = space->getNumKeys(ele);
 	int nbFF = ele->getNumShapeFunctions();
@@ -193,7 +193,7 @@ void dG3DNonLocalHomogenizedTangent::get(MElement *ele,int npts,IntPt *GP,fullMa
 void dG3DExtraDofHomogenizedTangent::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) const{
 	dG3DHomogenizedTangent::get(ele,npts,GP,m);
 	
-	FunctionSpaceBase* space = _dom->getFunctionSpace();
+	const FunctionSpaceBase* space = _dom->getFunctionSpace();
 	
   int nbDof = space->getNumKeys(ele);
   int nbFF = ele->getNumShapeFunctions();
diff --git a/dG3D/src/extraFieldFunctionSpace.h b/dG3D/src/extraFieldFunctionSpace.h
index 0aaeb2acceb79766cdab6894c3d74eebf1ae779e..9b4056c0a8df9b29b4398ec789ce71412cc9979b 100644
--- a/dG3D/src/extraFieldFunctionSpace.h
+++ b/dG3D/src/extraFieldFunctionSpace.h
@@ -17,7 +17,7 @@
 
 class extraFieldFunctionSpace : public ThreeDLagrangeFunctionSpace{
   protected:
-    virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys){
+    virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const{
       Msg::Fatal("This function is not implemented");
     }
 	
@@ -38,7 +38,7 @@ class extraFieldFunctionSpace : public ThreeDLagrangeFunctionSpace{
 			else if (_ncomp == 3) return new extraFieldFunctionSpace(_dom,id,3,comp[0],comp[1],comp[2]);
 			else return new extraFieldFunctionSpace(_dom,id,_ncomp);
 		}
-    virtual int getNumKeys(MElement *ele) {
+    virtual int getNumKeys(MElement *ele) const {
       if (_dom->IsMultipleFieldFormulation()){
 				if (dynamic_cast<MInterfaceElement*>(ele) == NULL) {
 					Msg::Fatal("Multiple fields are now considered for interface element only");
@@ -80,7 +80,7 @@ class extraFieldFunctionSpace : public ThreeDLagrangeFunctionSpace{
 			}
     };
 
-    virtual void getKeys(MElement *ele, std::vector<Dof> &keys) {
+    virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const{
 			if (_dom->IsMultipleFieldFormulation()){
 				if (_dom->getEnrichedUnknownLocation() == partDomain::NODE){
 					// generate Dof at nodes
diff --git a/dG3D/src/multipleFieldDG3DTerms.cpp b/dG3D/src/multipleFieldDG3DTerms.cpp
index f6f82f8b83bd42ac8a2bb91c103583f306fe3329..6fa6e1d88bb3df477b58001851c86ea8a0bd7747 100644
--- a/dG3D/src/multipleFieldDG3DTerms.cpp
+++ b/dG3D/src/multipleFieldDG3DTerms.cpp
@@ -15,7 +15,7 @@
 
 void twoFieldDG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFor) const
 {
-  FunctionSpaceBase* space  = _dom->getFunctionSpace();
+  const FunctionSpaceBase* space  = _dom->getFunctionSpace();
   int nbdof = space->getNumKeys(ele);
   int nbFF = ele->getNumShapeFunctions();
 	
@@ -59,7 +59,7 @@ void twoFieldDG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<doub
 
 void twoFieldDG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &mStiff) const
 {
-  FunctionSpaceBase* space  = _dom->getFunctionSpace();
+  const FunctionSpaceBase* space  = _dom->getFunctionSpace();
   int nbdof = space->getNumKeys(ele);
   int nbFF = ele->getNumShapeFunctions();
 	
@@ -143,9 +143,9 @@ void twoFieldDG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<
   const int nbFFp = ie->getElem(1)->getNumVertices();
   const int nbFF = ele->getNumVertices();
   
-  FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
-  FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
-  FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
+  const FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
+  const FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
+  const FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
 
   const int nbdofm = minusSpace->getNumKeys(ie->getElem(0));
   const int nbdofp = plusSpace->getNumKeys(ie->getElem(1));
@@ -296,9 +296,9 @@ void twoFieldDG3DStiffnessInter::get(MElement *ele,int npts,IntPt *GP,fullMatrix
   const int nbFFp = ie->getElem(1)->getNumVertices();
   const int nbFF = ele->getNumVertices();
 
-  FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
-  FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
-  FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
+  const FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
+  const FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
+  const FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
 
   const int nbdofm = minusSpace->getNumKeys(ie->getElem(0));
   const int nbdofp = plusSpace->getNumKeys(ie->getElem(1));
@@ -630,9 +630,9 @@ void dG3DTwoFieldForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<
   const int nbFFp = ie->getElem(1)->getNumVertices();
   const int nbFF = ele->getNumVertices();
   
-  FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
-  FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
-  FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
+  const FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
+  const FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
+  const FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
 
   const int nbdofm = minusSpace->getNumKeys(ie->getElem(0));
   const int nbdofp = plusSpace->getNumKeys(ie->getElem(1));
@@ -841,9 +841,9 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri
   const int nbFFp = ie->getElem(1)->getNumVertices();
   const int nbFF = ele->getNumVertices();
   
-  FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
-  FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
-  FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
+  const FunctionSpaceBase* minusSpace  = _dom->getFunctionSpaceMinus();
+  const FunctionSpaceBase* plusSpace = _dom->getFunctionSpacePlus();
+  const FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace(); 
 
   const int nbdofm = minusSpace->getNumKeys(ie->getElem(0));
   const int nbdofp = plusSpace->getNumKeys(ie->getElem(1));
diff --git a/dG3D/src/pathFollowingTerms.cpp b/dG3D/src/pathFollowingTerms.cpp
index 46b7e6f2314a31e4ff6f7a7f3f1ea9172267ab78..fff39643880498571dc431385fde476655e5379a 100644
--- a/dG3D/src/pathFollowingTerms.cpp
+++ b/dG3D/src/pathFollowingTerms.cpp
@@ -28,7 +28,7 @@ void dG3DDissipationPathFollowingBulkScalarTerm::get(MElement *ele, int npts, In
 };
 
 void dG3DDissipationPathFollowingBulkLinearTerm::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) const{
-	FunctionSpaceBase* sp = _dom->getFunctionSpace();
+	const FunctionSpaceBase* sp = _dom->getFunctionSpace();
 	
 	int nbFF = ele->getNumShapeFunctions();
 	int nbdof = sp->getNumKeys(ele);
@@ -67,7 +67,7 @@ void dG3DEnhancedStrainDissipationPathFollowingBulkLinearTerm::get(MElement *ele
 	
 	dG3DDissipationPathFollowingBulkLinearTerm::get(ele,npts,GP,m);
 	
-	FunctionSpaceBase* sp = _dom->getFunctionSpace();
+	const FunctionSpaceBase* sp = _dom->getFunctionSpace();
 	int nbFF = ele->getNumShapeFunctions();
 
 	const AllIPState::ipstateElementContainer *vips = _ipf->getAips()->getIPstate(ele->getNum());
@@ -140,9 +140,9 @@ void dG3DDissipationPathFollowingBoundScalarTerm::get(MElement *ele, int npts, I
 void dG3DDissipationPathFollowingBoundLinearTerm::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) const{
 	if(_dom->getFormulation())
 	{
-		FunctionSpaceBase* spaceMinus = _dom->getFunctionSpaceMinus();
-		FunctionSpaceBase* spacePlus = _dom->getFunctionSpacePlus();
-		FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace();
+		const FunctionSpaceBase* spaceMinus = _dom->getFunctionSpaceMinus();
+		const FunctionSpaceBase* spacePlus = _dom->getFunctionSpacePlus();
+		const FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace();
 		
 		MInterfaceElement *ie = dynamic_cast<MInterfaceElement*>(ele);
 		const int nbFFm = ie->getElem(0)->getNumVertices();
@@ -259,9 +259,9 @@ void twoFieldDG3DDissipationPathFollowingBoundLinearTerm::get(MElement *ele,int
 	{
 		MInterfaceElement *ie = dynamic_cast<MInterfaceElement*>(ele);
 		
-		FunctionSpaceBase* spaceMinus = _dom->getFunctionSpaceMinus();
-		FunctionSpaceBase* spacePlus = _dom->getFunctionSpacePlus();
-		FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace();
+		const FunctionSpaceBase* spaceMinus = _dom->getFunctionSpaceMinus();
+		const FunctionSpaceBase* spacePlus = _dom->getFunctionSpacePlus();
+		const FunctionSpaceBase* spaceInter = _dom->getInterfaceFunctionSpace();
 		
 		const int nbFFm = ie->getElem(0)->getNumVertices();
 		const int nbFFp = ie->getElem(1)->getNumVertices();
@@ -356,7 +356,7 @@ void dG3DNonLocalDissipationPathFollowingBulkLinearTerm::get(MElement *ele,int n
 	
 	const AllIPState::ipstateElementContainer *vips = _ipf->getAips()->getIPstate(ele->getNum());
 	
-	FunctionSpaceBase* space = _dom->getFunctionSpace();
+	const FunctionSpaceBase* space = _dom->getFunctionSpace();
 		
 	int nbFF = ele->getNumShapeFunctions();
 	for (int i = 0; i < npts; i++)
diff --git a/dgshell/src/dgC0ShellFunctionSpace.h b/dgshell/src/dgC0ShellFunctionSpace.h
index af45c19655c042a9e6db9104a068f39f761c63f3..4c6a95baefd043d5c147f377b4f49dcdd932b4fe 100644
--- a/dgshell/src/dgC0ShellFunctionSpace.h
+++ b/dgshell/src/dgC0ShellFunctionSpace.h
@@ -369,7 +369,7 @@ class DgC0FullDgLagrangeBetween2DomainsFunctionSpace : public IsoparametricLagra
  protected:
   FunctionSpaceBase *spaceMinus;
   FunctionSpaceBase *spacePlus;
-  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeysOnElement(MElement *ele, std::vector<Dof> &keys) const {
     Msg::Error("Impossible to get key for an element on a interface Domain");
   }
  public:
@@ -381,18 +381,21 @@ class DgC0FullDgLagrangeBetween2DomainsFunctionSpace : public IsoparametricLagra
 		return new DgC0FullDgLagrangeBetween2DomainsFunctionSpace(spaceMinus,spacePlus);
 	}
 
-  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys){
+  virtual void getKeys(MInterfaceElement *iele, std::vector<Dof> &keys) const{
     spaceMinus->getKeys(iele->getElem(0),keys);
     spacePlus->getKeys(iele->getElem(1),keys);
   }
-  virtual void getKeys(MElement *ele, std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele, std::vector<Dof> &keys) const{
     // As all element are interface element on an interface domain
     MInterfaceElement *iele = dynamic_cast<MInterfaceElement*>(ele);
     this->getKeys(iele,keys);
   }
   // special function of interFunctionSpace
-  virtual FunctionSpaceBase* getMinusSpace()const {return spaceMinus;}
-  virtual FunctionSpaceBase* getPlusSpace() const {return spacePlus;}
+  virtual FunctionSpaceBase* getMinusSpace() {return spaceMinus;}
+  virtual FunctionSpaceBase* getPlusSpace()  {return spacePlus;}
+  
+  virtual const FunctionSpaceBase* getMinusSpace()const {return spaceMinus;}
+  virtual const FunctionSpaceBase* getPlusSpace() const {return spacePlus;}
   virtual void getNumKeys(MInterfaceElement *ele, int &numMinus, int &numPlus) const
   {
     numMinus = spaceMinus->getNumKeys(ele->getElem(0));
diff --git a/dgshell/src/dgC0ShellRigidContactSpace.h b/dgshell/src/dgC0ShellRigidContactSpace.h
index db206eb398990065fcc7a5375d0fda5affdc9f21..d5df9836c2df0d2ceb0c3375db59469f4ccba467 100644
--- a/dgshell/src/dgC0ShellRigidContactSpace.h
+++ b/dgshell/src/dgC0ShellRigidContactSpace.h
@@ -18,7 +18,7 @@ class dgC0ShellRigidContactSpace : public rigidContactSpaceBase{
   dgC0ShellRigidContactSpace(const int id, FunctionSpace<double> *sp, MVertex *ver) : rigidContactSpaceBase(id,sp,ver){}
   dgC0ShellRigidContactSpace(const int id, FunctionSpace<double> *sp, MVertex *ver, const int comp1,
                    const int comp2 = -1, const int comp3 =-1) : rigidContactSpaceBase(id,sp,ver,comp1,comp2,comp3){}
-  virtual void getKeys(MElement *ele,std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele,std::vector<Dof> &keys) const{
     if(ele->getDim() != 0){ // if dim 0 return the key of gravity center only !!
       _space->getKeys(ele,keys);
     }
@@ -26,18 +26,18 @@ class dgC0ShellRigidContactSpace : public rigidContactSpaceBase{
       keys.push_back(Dof(_gc->getNum(),Dof3IntType::createTypeWithThreeInts(_comp[i],_id)));
     }
   }
-  virtual int getNumKeys(MElement *ele){
+  virtual int getNumKeys(MElement *ele) const{
     int nkeysele = _space->getNumKeys(ele);
     return nkeysele + _comp.size();
   }
-  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys){
+  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys) const{
     for(int i=0;i<_comp.size(); i++)
       keys.push_back(Dof(_gc->getNum(),Dof3IntType::createTypeWithThreeInts(_comp[i],_id)));
   }
-  virtual int getNumKeysOfGravityCenter(){
+  virtual int getNumKeysOfGravityCenter() const{
     return _comp.size();
   }
-  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys){
+  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys) const{
     // generate keys of element and select the good ones after LAGRANGE OK ?? CHANGE THIS HOW TODO
     std::vector<Dof> tkeys;
     this->getKeys(ele,tkeys);
@@ -48,7 +48,7 @@ class dgC0ShellRigidContactSpace : public rigidContactSpaceBase{
       keys.push_back(tkeys[ind+i*nbver]);
     this->getKeysOfGravityCenter(keys);
   }
-  virtual int getNumKeys(MElement *ele, int ind){
+  virtual int getNumKeys(MElement *ele, int ind) const{
     return 2*_comp.size();
   }
 };
diff --git a/dgshell/src/dgShellDomain.h b/dgshell/src/dgShellDomain.h
index 629951225279161516e88fd0c83cdf2eee3adb61..39592a6cf4e0b49183d0a38d1b35df4c7b7ab239 100644
--- a/dgshell/src/dgShellDomain.h
+++ b/dgshell/src/dgShellDomain.h
@@ -77,10 +77,17 @@ class dgLinearShellDomain : public dgPartDomain{
                                                     const simpleFunctionTime<double>* f, const unknownField *uf,
                                                     const IPField * ip, nonLinearBoundaryCondition::location onWhat) const;
   virtual double scaleTimeStep() const {return _sts;}
-  virtual FunctionSpaceBase* getFunctionSpace() const{return _space;}
-  virtual FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
-  virtual FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
-	virtual FunctionSpaceBase* getInterfaceFunctionSpace() const {
+  virtual FunctionSpaceBase* getFunctionSpace(){return _space;}
+  virtual const FunctionSpaceBase* getFunctionSpace() const{return _space;}
+  virtual FunctionSpaceBase* getFunctionSpaceMinus(){return _spaceMinus;}
+  virtual const FunctionSpaceBase* getFunctionSpaceMinus() const{return _spaceMinus;}
+  virtual FunctionSpaceBase* getFunctionSpacePlus() {return _spacePlus;}
+  virtual const FunctionSpaceBase* getFunctionSpacePlus() const{return _spacePlus;}
+	virtual FunctionSpaceBase* getInterfaceFunctionSpace() {
+		Msg::Error("getInterfaceFunctionSpace is not defined");
+		return NULL;
+	}
+  virtual const FunctionSpaceBase* getInterfaceFunctionSpace() const {
 		Msg::Error("getInterfaceFunctionSpace is not defined");
 		return NULL;
 	}