diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 38e33dacd1de6e2d48cf6c9bd0c35bc8ddb42133..2d2599c0b5575956f6302d6172fab4efcdedb261 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -429,21 +429,21 @@ void MyelasticitySolver::solve()
 
   for (unsigned int i = 0; i < edgeNeu.size(); i++)
   {
-    LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,edgeNeu[i]._f);
+    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);
   }
 
   for (unsigned int i = 0; i < faceNeu.size(); i++)
   {
-    LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,faceNeu[i]._f);
+    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);
   }
 
   for (unsigned int i = 0; i < volumeNeu.size(); i++)
   {
-    LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,volumeNeu[i]._f);
+    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);
   }
@@ -451,7 +451,7 @@ void MyelasticitySolver::solve()
   GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
   for (unsigned int i = 0; i < elasticFields.size(); i++)
   {
-    ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
+    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);
   }
 
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index a86e28117ae0bb184fb6348477f58d02188b2fe1..70c1abb84860550eb746eeb4fe194641aad9efb5 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -19,7 +19,7 @@
 #include "terms.h"
 #include "quadratureRules.h"
 #include "MVertex.h"
-#include <iostream>
+
 
 template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler) // symmetric
 {
diff --git a/Solver/terms.h b/Solver/terms.h
index 6ac36ff7d6de9807e50d1af3e09177580162c411..721574d943ab6c014b66654e3129503a51cce76a 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -105,7 +105,7 @@ class ScalarTerm : public ScalarTermBase
 
 
 
-template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // non symmetric 
+template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,S2> // non symmetric 
 {
  protected : 
   double E,nu;
@@ -118,7 +118,7 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no
       {  0,   0,   0,    0,   0, C44} };*/
 
  public : 
-  ElasticTerm(S1& space1_,S2& space2_,double E_,double nu_) : BilinearTerm<S1,S2>(space1_,space2_),E(E_),nu(nu_),H(6,6)
+  IsotropicElasticTerm(S1& space1_,S2& space2_,double E_,double nu_) : BilinearTerm<S1,S2>(space1_,space2_),E(E_),nu(nu_),H(6,6)
   {
     double FACT = E / (1 + nu);
     double C11 = FACT * (1 - nu) / (1 - 2 * nu);
@@ -176,7 +176,7 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no
 
 
 
-template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
+template<class S1> class IsotropicElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
 {
  protected : 
   double E,nu;
@@ -189,7 +189,7 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
       {  0,   0,   0,    0,   0, C44} };*/
 
  public : 
-  ElasticTerm(S1& space1_,double E_,double nu_) : BilinearTerm<S1,S1>(space1_,space1_),E(E_),nu(nu_),H(6,6)
+  IsotropicElasticTerm(S1& space1_,double E_,double nu_) : BilinearTerm<S1,S1>(space1_,space1_),E(E_),nu(nu_),H(6,6)
   {
     double FACT = E / (1 + nu);
     double C11 = FACT * (1 - nu) / (1 - 2 * nu);
@@ -235,36 +235,12 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
 inline double dot(const double &a, const double &b)
 { return a*b; }
 
-template<class S1> class LoadTerm : public LinearTerm<S1>
-{
-  typename S1::ValType Load;
- public : 
-  LoadTerm(S1& space1_,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);
-    double jac[3][3];
-    m.resize(nbFF);
-    m.scale(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::ValType> Vals;
-      LinearTerm<S1>::space1.f(ele,u, v, w, Vals);
-      for (int j = 0; j < nbFF ; ++j)
-      {
-        m(j)+=dot(Vals[j],Load)*weight*detJ;
-      }
-    }
-  }
-};
 
-template<class S1> class LoadTermSF : public LinearTerm<S1>
+template<class S1> class LoadTerm : public LinearTerm<S1>
 {
   simpleFunction<typename S1::ValType> Load;
  public : 
-  LoadTermSF(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);