less template, more inheritance

Solver/elasticitySolver.cpp

@@ -188,7 +188,7 @@ void elasticitySolver::solve()
 
   for (unsigned int i = 0; i < allNeumann.size(); i++)
   {
-    LoadTerm<FunctionSpace<SVector3> > Lterm(*LagSpace,allNeumann[i]._f);
+    LoadTerm<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
@@ -200,13 +200,23 @@ void elasticitySolver::solve()
   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);
+    IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
+//    LaplaceTerm<SVector3,SVector3> Eterm(*LagSpace);
     Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
   }
 
   printf("-- done assembling!\n");
   lsys->systemSolve();
   printf("-- done solving!\n");
+  double energ=0;
+  for (unsigned int i = 0; i < elasticFields.size(); i++)
+  {
+    SolverField<SVector3> Field(pAssembler, LagSpace);
+    IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
+  }
+  printf("elastic energy=%f\n",energ);
 }
 
 
@@ -264,21 +274,12 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
       for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
     }
   }
-  std::map<int, std::vector<double> > data;  
+  std::map<int, std::vector<double> > data;
+  SolverField<SVector3> Field(pAssembler, LagSpace);
   for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
   {
-    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];
-    }
+    SVector3 val;
+    Field.f(*it,val);
     std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
     data[(*it)->getNum()]=vec;
   }
@@ -290,30 +291,54 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
 {
   std::map<int, std::vector<double> > data;
   GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
-  double energ=0;
   for (unsigned int i = 0; i < elasticFields.size(); ++i)
   {
     SolverField<SVector3> Field(pAssembler, LagSpace);
-    IsotropicElasticTerm<SolverField<SVector3> ,SolverField<SVector3> > Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
-    ScalarTermOne One;
+    IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    ScalarTermConstant One(1.0);
     for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
     {
       MElement *e=*it;
-      fullMatrix<double> localMatrix;
+      double energ;
+      double vol;
       IntPt *GP;
       int npts=Integ_Bulk.getIntPoints(e,&GP);
-      Eterm.get(e,npts,GP,localMatrix);
-      double vol=0;
+      Elastic_Energy_Term.get(e,npts,GP,energ);
       One.get(e,npts,GP,vol);
       std::vector<double> vec;
-      vec.push_back(localMatrix(0,0)/vol);
-      energ+=localMatrix(0,0);
+      vec.push_back(energ/vol);
       data[(*it)->getNum()]=vec;
     }
   }
-  std::cout<< "elastic energy=" << energ << std::endl;
   PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
-  return pv;  
+  return pv;
+}
+
+PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
+{
+  std::map<int, std::vector<double> > data;
+  GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+  for (unsigned int i = 0; i < elasticFields.size(); ++i)
+  {
+    SolverField<SVector3> Field(pAssembler, LagSpace);
+    IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
+    {
+      MElement *e=*it;
+      double energ;
+      double vol;
+      IntPt *GP;
+      int npts=Integ_Bulk.getIntPoints(e,&GP);
+      Elastic_Energy_Term.get(e,npts,GP,energ);
+      std::vector<double> vec;
+      vec.push_back(energ);
+      data[(*it)->getNum()]=vec;
+    }
+  }
+  PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
+  return pv;
 }
 
 

Solver/elasticitySolver.h

@@ -20,9 +20,8 @@ class groupOfElements;
 struct elasticField {
   int _tag; // tag for the dofManager
   groupOfElements *g; // support for this field
-  simpleFunction<double> *_enrichment; // XFEM enrichment
   double _E, _nu; // specific elastic datas (should be somewhere else)
-  elasticField () : g(0), _enrichment(0),_tag(0){}
+  elasticField () : g(0),_tag(0){}
 };
 
 struct BoundaryCondition
@@ -75,7 +74,7 @@ class elasticitySolver
   virtual void solve();
   virtual PView *buildDisplacementView(const std::string &postFileName);
   virtual PView *buildElasticEnergyView(const std::string &postFileName);
-  // virtual 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

Solver/materialLaw.h

+//
+// C++ Interface: materialLaw
+//
+// Description: 
+//
+//
+// Author:  <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _MATERIALLAW_H_
+#define _MATERIALLAW_H_
+
+class Material
+{
+ public:
+  virtual ~Material() {}
+  
+};
+
+
+
+#endif //_MATERIALLAW_H_

Solver/solverAlgorithms.h

@@ -90,6 +90,30 @@ template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &
   assembler.assemble(R, localVector);
 }
 
+template<class Iterator,class dataMat> void Assemble(ScalarTermBase &term,Iterator itbegin,Iterator itend,QuadratureBase &integrator,dataMat & val)
+{
+  dataMat localval;
+  for (Iterator it = itbegin;it!=itend; ++it)
+  {
+    MElement *e = *it;
+    IntPt *GP;
+    int npts=integrator.getIntPoints(e,&GP);
+    term.get(e,npts,GP,localval);
+    val+=localval;
+  }
+}
+
+template<class Iterator,class dataMat> void Assemble(ScalarTermBase &term,MElement *e,QuadratureBase &integrator,dataMat & val)
+{
+  dataMat localval;
+  IntPt *GP;
+  int npts=integrator.getIntPoints(e,&GP);
+  term.get(e,npts,GP,localval);
+  val+=localval;
+}
+
+
+
 template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &dofs,std::vector<typename Assembler::dataVec> &vals)
 {
   int nbff=dofs.size();

Solver/terms.h

@@ -20,51 +20,77 @@
 #include "Numeric.h"
 #include "functionSpace.h"
 #include "groupOfElements.h"
+#include "materialLaw.h"
 
 class  BilinearTermBase
 {
  public :
+  virtual ~BilinearTermBase() {}
   virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
 };
 
-template<class S1,class S2> class BilinearTerm : public BilinearTermBase
+template<class T1,class T2> class BilinearTerm : public BilinearTermBase
 {
  protected :
-  S1& space1;
-  S2& space2;
+  FunctionSpace<T1>& space1;
+  FunctionSpace<T2>& space2;
  public :
-  BilinearTerm(S1& space1_,S2& space2_) : space1(space1_),space2(space2_) {}
+  BilinearTerm(FunctionSpace<T1>& space1_,FunctionSpace<T1>& space2_) : space1(space1_),space2(space2_) {}
+  virtual ~BilinearTerm() {}
 };
 
 class  LinearTermBase
 {
   public:
+  virtual ~LinearTermBase() {}
   virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
   virtual void get(MVertex *ver,fullVector<double> &m) =0;
 };
 
-template<class S1> class LinearTerm : public LinearTermBase
+template<class T1> class LinearTerm : public LinearTermBase
 {
  protected :
-  S1& space1;
- public : 
-  LinearTerm(S1& space1_) : space1(space1_) {}
+  FunctionSpace<T1>& space1;
+ public :
+  LinearTerm(FunctionSpace<T1>& space1_) : space1(space1_) {}
+  virtual ~LinearTerm() {}
 };
 
 class  ScalarTermBase
 {
  public : 
+  virtual ~ScalarTermBase() {}
   virtual void get(MElement *ele,int npts,IntPt *GP,double &val) =0;
 };
 
 class ScalarTerm : public ScalarTermBase
 {
- public : 
+ public :
+  virtual ~ScalarTerm() {}
 };
 
-class ScalarTermOne : public ScalarTerm
+
+template<class T1,class T2> class BilinearTermToScalarTerm : public ScalarTerm
 {
- public : 
+  BilinearTerm<T1,T2> &bilterm;
+  public :
+  BilinearTermToScalarTerm(BilinearTerm<T1,T2> &bilterm_): bilterm(bilterm_){}
+  virtual ~BilinearTermToScalarTerm() {}
+  virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
+  {
+    fullMatrix<double> localMatrix;
+    bilterm.get(ele,npts,GP,localMatrix);
+    val=localMatrix(0,0);
+  }
+};
+
+
+class ScalarTermConstant : public ScalarTerm
+{
+  double val;
+ public :
+  ScalarTermConstant(double val_=1.0): val(val_) {}
+  virtual ~ScalarTermConstant() {}
   virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
   {
     double jac[3][3];
@@ -76,16 +102,20 @@ class ScalarTermOne : public ScalarTerm
       val+=weight*detJ;
     }
   }
+  virtual void get(MVertex *ver,double &val)
+  {
+      val=1;
+  }
 };
 
 
 
-template<class S1,class S2> class LaplaceTerm : public BilinearTerm<S1,S2> 
+template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2> 
 {
  public : 
-  LaplaceTerm(S1& space1_,S2& space2_) : BilinearTerm<S1,S2>(space1_,space2_)
+  LaplaceTerm(FunctionSpace<T1>& space1_,FunctionSpace<T2>& space2_) : BilinearTerm<T1,T2>(space1_,space2_)
   {}
-  
+  virtual ~LaplaceTerm() {}
   virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
   {
     Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
@@ -98,15 +128,15 @@ template<class S1,class S2> class LaplaceTerm : public BilinearTerm<S1,S2>
 
 
 
-template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
+template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symmetric
 {
  public : 
-  LaplaceTerm(S1& space1_) : BilinearTerm<S1,S1>(space1_,space1_)
+  LaplaceTerm(FunctionSpace<T1>& space1_) : BilinearTerm<T1,T1>(space1_,space1_)
   {}
-  
+  virtual ~LaplaceTerm() {}
   virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
   {
-    int nbFF = BilinearTerm<S1,S1>::space1.getNumKeys(ele);
+    int nbFF = BilinearTerm<T1,T1>::space1.getNumKeys(ele);
     double jac[3][3];
     m.resize(nbFF, nbFF);
     m.setAll(0.);
@@ -114,8 +144,8 @@ template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
     {
       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);
+      std::vector<typename TensorialTraits<T1>::GradType> Grads;
+      BilinearTerm<T1,T1>::space1.gradf(ele,u, v, w, Grads);
       for (int j = 0; j < nbFF; j++)
       {
         for (int k = j; k < nbFF; k++)
@@ -134,10 +164,11 @@ template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
 
 
 
-template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,S2> // non symmetric 
+class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3> 
 {
  protected : 
   double E,nu;
+  bool sym;
   fullMatrix<double> H;/* =
     { {C11, C12, C12,    0,   0,   0},
       {C12, C11, C12,    0,   0,   0},
@@ -147,7 +178,7 @@ template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,
       {  0,   0,   0,    0,   0, C44} };*/
 
  public : 
-  IsotropicElasticTerm(S1& space1_,S2& space2_,double E_,double nu_) : BilinearTerm<S1,S2>(space1_,space2_),E(E_),nu(nu_),H(6,6)
+  IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,FunctionSpace<SVector3>& space2_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space2_),E(E_),nu(nu_),H(6,6)
   {
     double FACT = E / (1 + nu);
     double C11 = FACT * (1 - nu) / (1 - 2 * nu);
@@ -156,69 +187,9 @@ template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,
     H.scale(0.);
     for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
     H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
+    sym=(&space1_==&space2_);
   }
-  virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
-  {
-    int nbFF1 = BilinearTerm<S1,S2>::space1.getNumKeys(ele);
-    int nbFF2 = BilinearTerm<S1,S2>::space2.getNumKeys(ele);
-    double jac[3][3];
-    fullMatrix<double> B(6, nbFF2);
-    fullMatrix<double> BTH(nbFF2, 6);
-    fullMatrix<double> BT(nbFF1, 6);
-    m.resize(nbFF1, nbFF2);
-    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;// tableau de matrices...
-      std::vector<typename S2::GradType> GradsT;// tableau de matrices...
-      BilinearTerm<S1,S2>::space1.gradf(ele,u, v, w, Grads);
-      BilinearTerm<S1,S2>::space2.gradf(ele,u, v, w, GradsT);
-      for (int j = 0; j < nbFF1; j++)
-      {
-        BT(j, 0) = Grads[j](0,0);
-        BT(j, 1) = Grads[j](1,1);
-        BT(j, 2) = Grads[j](2,2);
-        BT(j, 3) = Grads[j](0,1)+Grads[j](1,0);
-        BT(j, 4) = Grads[j](1,2)+Grads[j](2,1);
-        BT(j, 5) = Grads[j](0,2)+Grads[j](2,0);
-      }
-      for (int j = 0; j < nbFF2; j++)
-      {
-        B(0, j) = GradsT[j](0,0);
-        B(1, j) = GradsT[j](1,1);
-        B(2, j) = GradsT[j](2,2);
-        B(3, j) = GradsT[j](0,1)+GradsT[j](1,0);
-        B(4, j) = GradsT[j](1,2)+GradsT[j](2,1);
-        B(5, j) = GradsT[j](0,2)+GradsT[j](2,0);
-      }
-      BTH.setAll(0.);
-      BTH.gemm(BT, H);
-      m.gemm(BTH, B, weight * detJ, 1.);
-    }
-  }
-}; // class
-
-
-
-
-
-
-template<class S1> class IsotropicElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
-{
- protected : 
-  double E,nu;
-  fullMatrix<double> H;/* =
-    { {C11, C12, C12,    0,   0,   0},
-      {C12, C11, C12,    0,   0,   0},
-      {C12, C12, C11,    0,   0,   0},
-      {  0,   0,   0,  C44,   0,   0},
-      {  0,   0,   0,    0, C44,   0},
-      {  0,   0,   0,    0,   0, C44} };*/
-
- public : 
-  IsotropicElasticTerm(S1& space1_,double E_,double nu_) : BilinearTerm<S1,S1>(space1_,space1_),E(E_),nu(nu_),H(6,6)
+  IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space1_),E(E_),nu(nu_),H(6,6)
   {
     double FACT = E / (1 + nu);
     double C11 = FACT * (1 - nu) / (1 - 2 * nu);
@@ -227,35 +198,80 @@ template<class S1> class IsotropicElasticTerm<S1,S1> : public BilinearTerm<S1,S1
     H.scale(0.);
     for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
     H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
+    sym=true;
   }
-  
+  virtual ~IsotropicElasticTerm() {}
   virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
   {
-    int nbFF = BilinearTerm<S1,S1>::space1.getNumKeys(ele);
-    double jac[3][3];
-    fullMatrix<double> B(6, nbFF);
-    fullMatrix<double> BTH(nbFF, 6);
-    fullMatrix<double> BT(nbFF, 6);
-    m.resize(nbFF, nbFF);
-    m.setAll(0.);
-    for (int i = 0; i < npts; i++)
+    if (sym)
     {
-      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); // a optimiser ??
-      for (int j = 0; j < nbFF; j++)
+      int nbFF = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
+      double jac[3][3];
+      fullMatrix<double> B(6, nbFF);
+      fullMatrix<double> BTH(nbFF, 6);
+      fullMatrix<double> BT(nbFF, 6);
+      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<TensorialTraits<SVector3>::GradType> Grads;
+        BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads); // a optimiser ??
+        for (int j = 0; j < nbFF; j++)
+        {
+          BT(j, 0) = B(0, j) = Grads[j](0,0);
+          BT(j, 1) = B(1, j) = Grads[j](1,1);
+          BT(j, 2) = B(2, j) = Grads[j](2,2);
+          BT(j, 3) = B(3, j) = Grads[j](0,1)+Grads[j](1,0);
+          BT(j, 4) = B(4, j) = Grads[j](1,2)+Grads[j](2,1);
+          BT(j, 5) = B(5, j) = Grads[j](0,2)+Grads[j](2,0);
+        }
+        BTH.setAll(0.);
+        BTH.gemm(BT, H);
+        m.gemm(BTH, B, weight * detJ, 1.);
+      }
+    }
+    else
+    {
+      int nbFF1 = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
+      int nbFF2 = BilinearTerm<SVector3,SVector3>::space2.getNumKeys(ele);
+      double jac[3][3];
+      fullMatrix<double> B(6, nbFF2);
+      fullMatrix<double> BTH(nbFF2, 6);
+      fullMatrix<double> BT(nbFF1, 6);
+      m.resize(nbFF1, nbFF2);
+      m.setAll(0.);
+      for (int i = 0; i < npts; i++)
       {
-        BT(j, 0) = B(0, j) = Grads[j](0,0);
-        BT(j, 1) = B(1, j) = Grads[j](1,1);
-        BT(j, 2) = B(2, j) = Grads[j](2,2);
-        BT(j, 3) = B(3, j) = Grads[j](0,1)+Grads[j](1,0);
-        BT(j, 4) = B(4, j) = Grads[j](1,2)+Grads[j](2,1);
-        BT(j, 5) = B(5, j) = Grads[j](0,2)+Grads[j](2,0);
+        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<TensorialTraits<SVector3>::GradType> Grads;// tableau de matrices...
+        std::vector<TensorialTraits<SVector3>::GradType> GradsT;// tableau de matrices...
+        BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads);
+        BilinearTerm<SVector3,SVector3>::space2.gradf(ele,u, v, w, GradsT);
+        for (int j = 0; j < nbFF1; j++)
+        {
+          BT(j, 0) = Grads[j](0,0);
+          BT(j, 1) = Grads[j](1,1);
+          BT(j, 2) = Grads[j](2,2);
+          BT(j, 3) = Grads[j](0,1)+Grads[j](1,0);
+          BT(j, 4) = Grads[j](1,2)+Grads[j](2,1);
+          BT(j, 5) = Grads[j](0,2)+Grads[j](2,0);
+        }
+        for (int j = 0; j < nbFF2; j++)
+        {
+          B(0, j) = GradsT[j](0,0);
+          B(1, j) = GradsT[j](1,1);
+          B(2, j) = GradsT[j](2,2);
+          B(3, j) = GradsT[j](0,1)+GradsT[j](1,0);
+          B(4, j) = GradsT[j](1,2)+GradsT[j](2,1);
+          B(5, j) = GradsT[j](0,2)+GradsT[j](2,0);
+        }
+        BTH.setAll(0.);
+        BTH.gemm(BT, H);
+        m.gemm(BTH, B, weight * detJ, 1.);
       }
-      BTH.setAll(0.);
-      BTH.gemm(BT, H);
-      m.gemm(BTH, B, weight * detJ, 1.);
     }
   }
 }; // class
@@ -265,14 +281,16 @@ inline double dot(const double &a, const double &b)
 { return a*b; }
 
 
-template<class S1> class LoadTerm : public LinearTerm<S1>
+template<class T1> class LoadTerm : public LinearTerm<T1>
 {
-  simpleFunction<typename S1::ValType> &Load;
+  simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
  public : 
-  LoadTerm(S1& space1_,simpleFunction<typename S1::ValType> &Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
+  LoadTerm(FunctionSpace<T1>& space1_,simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :LinearTerm<T1>(space1_),Load(Load_) {}
+  virtual ~LoadTerm() {}
+
   virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
   {
-    double nbFF=LinearTerm<S1>::space1.getNumKeys(ele);
+    double nbFF=LinearTerm<T1>::space1.getNumKeys(ele);
     double jac[3][3];
     m.resize(nbFF);
     m.scale(0.);
@@ -280,11 +298,11 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
     {
       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);
+      std::vector<typename TensorialTraits<T1>::ValType> Vals;
+      LinearTerm<T1>::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());
+      typename TensorialTraits<T1>::ValType load=Load(p.x(),p.y(),p.z());
       for (int j = 0; j < nbFF ; ++j)
       {
         m(j)+=dot(Vals[j],load)*weight*detJ;
@@ -294,12 +312,12 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
 
   virtual void get(MVertex *ver,fullVector<double> &m)
   {
-    double nbFF=LinearTerm<S1>::space1.getNumKeys(ver);
+    double nbFF=LinearTerm<T1>::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());
+    std::vector<typename TensorialTraits<T1>::ValType> Vals;
+    LinearTerm<T1>::space1.f(ver, Vals);
+    typename TensorialTraits<T1>::ValType load=Load(ver->x(),ver->y(),ver->z());
     for (int j = 0; j < nbFF ; ++j)
     {
       m(j)=dot(Vals[j],load);