diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 21b20c8590d4c8339524bdbbc818d92bf6c21c2d..117005fdcfd439abaeaac523e5de690773dd17a9 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -513,38 +513,9 @@ void MyelasticitySolver::solve()
   for (unsigned int i = 0; i < elasticFields.size(); i++)
   {
     ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
-    fullMatrix<double> localMatrix;
-    std::vector<Dof> R;
-    for ( groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end() ; ++it)
-    {
-      MElement *e = *it;
-      R.clear();
-      int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
-      int npts=0;
-      IntPt *GP;
-      e->getIntegrationPoints(integrationOrder, &npts, &GP);
-      Eterm.get(e,npts,GP,localMatrix);
-      P123.getKeys(e,R);
-      pAssembler->assemble(R, localMatrix);
-    }
+    Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),*pAssembler);
   }
 
-
-/*
-
-  for (std::map<int, SVector3 >::iterator it = volumeForces.begin();it != volumeForces.end(); ++it)
-  {
-    int iVolume = it->first;
-    SVector3 f = it->second;
-    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]);
-    }
-  }
-*/
 /*
     for (int i=0;i<pAssembler->sizeOfR();i++)
     {
diff --git a/Solver/terms.h b/Solver/terms.h
index b229a2f5556c15449b7504035505228e89092165..62161ba0612ad4f7972143b4411d9640396268b8 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -259,5 +259,22 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
 };
 
 
+template<class T,class S1, class I,class A> void Assemble(T& term,S1& space1,I itbegin,I itend,A& assembler) // symmetric
+{
+    fullMatrix<double> localMatrix;
+    std::vector<Dof> R;
+    for (I it = itbegin;it!=itend; ++it)
+    {
+      MElement *e = *it;
+      R.clear();
+      int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
+      int npts=0;
+      IntPt *GP;
+      e->getIntegrationPoints(integrationOrder, &npts, &GP);
+      term.get(e,npts,GP,localMatrix);
+      space1.getKeys(e,R);
+      assembler.assemble(R, localMatrix);
+    }
+}
 
 #endif// _TERMS_H_