No commit message

No commit message

Solver/CMakeLists.txt

@@ -19,6 +19,8 @@ set(SRC
   functionSpace.cpp
   filters.cpp
   sparsityPattern.cpp
+  STensor43.cpp
+  terms.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) 

Solver/STensor43.cpp

+
+#include "STensor43.h"
+
+void STensor43::print (const char *s) const
+{
+  printf("STensor43::print to be implemented \n");
+/*  printf(" tensor4 %s : \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n",s,
+         (*this)(0,0),(*this)(0,1),(*this)(0,2),
+         (*this)(1,0),(*this)(1,1),(*this)(1,2),
+         (*this)(2,0),(*this)(2,1),(*this)(2,2));*/
+}

Solver/STensor43.h

+#ifndef _STENSOR43_H_
+#define _STENSOR43_H_
+
+#include "STensor3.h"
+#include "fullMatrix.h"
+#include "Numeric.h"
+
+
+// concrete class for general 3x3 matrix
+
+class STensor43 {
+ protected:
+  // 0000 0001 0002 0010 ... 2211 2212 2220 2221 2222 
+  double _val[81];
+ public:
+  inline int getIndex(int i, int j, int k, int l) const
+  {
+    static int _index[3][3][3][3] = {{{{0,1,2},{3,4,5},{6,7,8}},{{9,10,11},{12,13,14},{15,16,17}},{{18,19,20},{21,22,23},{24,25,26}}},
+                                     {{{27,28,29},{30,31,32},{33,34,35}},{{36,37,38},{39,40,41},{42,43,44}},{{45,46,47},{48,49,50},{51,52,53}}},
+                                     {{{54,55,56},{57,58,59},{60,61,62}},{{63,64,65},{66,67,68},{69,70,71}},{{72,73,74},{75,76,77},{77,79,80}}}};
+    return _index[i][j][k][l];
+  }
+  STensor43(const STensor43& other)
+  {
+    for (int i = 0; i < 81; i++) _val[i] = other._val[i];
+  }
+  // default constructor, null tensor
+  STensor43(const double v = 0.0)
+  {
+    for (int i = 0; i < 3; i++)
+      for (int j = 0; j < 3; j++)
+        for (int k = 0; k < 3; k++)
+          for (int l = 0; l < 3; l++)
+            if ((i==k)&&(j==l))
+              _val[getIndex(i, j, k, l)]=v;
+            else 
+              _val[getIndex(i, j, k, l)]=0.0;
+  }
+  inline double &operator()(int i, int j,int k, int l)
+  {
+    return _val[getIndex(i, j, k, l)];
+  }
+  inline double operator()(int i, int j, int k, int l) const
+  {
+    return _val[getIndex(i, j, k ,l)];
+  }  
+  STensor43 operator + (const STensor43 &other) const
+  {
+    STensor43 res(*this);
+    for (int i = 0; i < 81; i++) res._val[i] += other._val[i];
+    return res;
+  }
+  STensor43& operator += (const STensor43 &other)
+  {
+    for (int i = 0; i < 81; i++) _val[i] += other._val[i];
+    return *this;
+  }
+  STensor43& operator *= (const double &other) 
+  {
+    for (int i = 0; i < 81; i++) _val[i] *= other;
+    return *this;
+  }
+/*  STensor43& operator *= (const STensor43 &other) 
+  {
+// to be implemented
+    return *this;
+  }*/
+  void print(const char *) const;
+};
+
+// tensor product
+inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c)
+{ 
+    for (int i = 0; i < 3; i++)
+      for (int j = 0; j < 3; j++)
+        for (int k = 0; k < 3; k++)
+          for (int l = 0; l < 3; l++)
+            c(i,j,k,l)=a(i,j)*b(k,l);
+}
+
+inline double dot(const STensor43 &a, const STensor43 &b)
+{ 
+  double prod=0;  
+  for (int i = 0; i < 3; i++)
+    for (int j = 0; j < 3; j++)
+      for (int k = 0; k < 3; k++)
+        for (int l = 0; l < 3; l++)
+      prod+=a(i,j,k,l)*b(i,j,k,l);
+  return prod;
+}
+
+inline STensor43 operator*(const STensor43 &t, double m)
+{ 
+  STensor43 val(t);
+  val *= m;
+  return val;
+}
+
+inline STensor43 operator*(double m,const STensor43 &t)
+{ 
+  STensor43 val(t);
+  val *= m;
+  return val;
+}
+
+inline STensor3 operator*(const STensor43 &t, const STensor3 &m)
+{ 
+  STensor3 val(0.);
+  for (int i = 0; i < 3; i++)
+    for (int j = 0; j < 3; j++)
+      for (int k = 0; k < 3; k++)
+        for (int l = 0; l < 3; l++)
+          val(i,j)+=t(i,j,k,l)*m(k,l);
+  return val;
+}
+
+inline STensor3 operator*( const STensor3 &m , const STensor43 &t)
+{ 
+  STensor3 val(0.);
+  for (int i = 0; i < 3; i++)
+    for (int j = 0; j < 3; j++)
+      for (int k = 0; k < 3; k++)
+        for (int l = 0; l < 3; l++)
+          val(k,l)+=t(i,j,k,l)*m(i,j);
+  return val;
+}
+
+
+
+#endif

Solver/elasticitySolver.cpp

@@ -92,7 +92,7 @@ void elasticitySolver::solve()
   {
     SolverField<SVector3> Field(pAssembler, LagSpace);
     IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
-    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    BilinearTermToScalarTerm 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);
@@ -107,7 +107,7 @@ void elasticitySolver::postSolve()
   {
     SolverField<SVector3> Field(pAssembler, LagSpace);
     IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
-    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    BilinearTermToScalarTerm 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);
@@ -551,7 +551,7 @@ PView* elasticitySolver::buildStressesView (const std::string postFileName)
     double nu = elasticFields[i]._nu;
     SolverField<SVector3> Field(pAssembler, LagSpace);
     IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
-    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
     for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
     {
       MElement *e=*it;
@@ -644,8 +644,8 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string postFileName)
     if(elasticFields[i]._E == 0.) continue;
     SolverField<SVector3> Field(pAssembler, LagSpace);
     IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
-    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
-    ScalarTermConstant One(1.0);
+    BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
+    ScalarTermConstant<double> One(1.0);
     for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
     {
       MElement *e = *it;
@@ -673,7 +673,7 @@ PView *elasticitySolver::buildVonMisesView(const std::string postFileName)
   {
     SolverField<SVector3> Field(pAssembler, LagSpace);
     IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
-    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
     for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
     {
       MElement *e=*it;

Solver/functionSpace.h

@@ -3,6 +3,7 @@
 
 #include "SVector3.h"
 #include "STensor3.h"
+#include "STensor43.h"
 #include <vector>
 #include <iterator>
 #include <iostream>
@@ -16,8 +17,8 @@ template<class T> struct TensorialTraits
 {
   typedef T ValType;
   typedef T GradType[3];
-/*  typedef T HessType[3][3];
-  typedef SVoid DivType;
+  typedef T HessType[3][3];
+/*  typedef SVoid DivType;
   typedef SVoid CurlType;*/
 };
 
@@ -26,6 +27,7 @@ template<> struct TensorialTraits<double>
   typedef double ValType;
   typedef SVector3 GradType;
   typedef STensor3 HessType;
+  typedef double TensProdType;
 /*  typedef SVoid DivType;
   typedef SVoid CurlType;*/
 };
@@ -35,10 +37,23 @@ template<> struct TensorialTraits<SVector3>
   typedef SVector3 ValType;
   typedef STensor3 GradType;
   typedef STensor3 HessType;
-  typedef double DivType;
-  typedef SVector3 CurlType;
+  typedef STensor3 TensProdType;
+//  typedef double DivType;
+//  typedef SVector3 CurlType;
 };
 
+template<> struct TensorialTraits<STensor3>
+{
+  typedef STensor3 ValType;
+//  typedef STensor3 GradType;
+//  typedef STensor3 HessType;
+//  typedef STensor3 TensProdType;
+  typedef STensor43 TensProdType;
+//  typedef double DivType;
+//  typedef SVector3 CurlType;
+};
+
+
 class FunctionSpaceBase
 {
  public:

Solver/solverAlgorithms.h

@@ -97,7 +97,7 @@ template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term,
 
 
 
-template<class Iterator, class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase &space,
+template<class Iterator, class Assembler> void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space,
                                                         Iterator itbegin, Iterator itend,
                                                         QuadratureBase &integrator, Assembler &assembler)
 {
@@ -114,7 +114,7 @@ template<class Iterator, class Assembler> void Assemble(LinearTermBase &term, Fu
   }
 }
 
-template<class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase &space, MElement *e,
+template<class Assembler> void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space, MElement *e,
                                         QuadratureBase &integrator, Assembler &assembler)
 {
   fullVector<typename Assembler::dataMat> localVector;
@@ -126,7 +126,7 @@ template<class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase
   assembler.assemble(R, localVector);
 }
 
-template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term,
+template<class Iterator, class dataMat> void Assemble(ScalarTermBase<double> &term,
                                                       Iterator itbegin, Iterator itend,
                                                       QuadratureBase &integrator, dataMat & val)
 {
@@ -140,7 +140,7 @@ template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term,
   }
 }
 
-template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term, MElement *e,
+template<class Iterator, class dataMat> void Assemble(ScalarTermBase<double> &term, MElement *e,
                                                       QuadratureBase &integrator, dataMat & val)
 {
   dataMat localval;

Solver/terms.cpp

+//
+// C++ Implementation: terms
+//
+// Description:
+//
+//
+// Author:  <Eric Bechet>, (C) 2011
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "terms.h"
+
+
+void BilinearTermToScalarTerm::get(MElement *ele, int npts, IntPt *GP, double &val) const
+{
+  fullMatrix<double> localMatrix;
+  bilterm.get(ele, npts, GP, localMatrix);
+  val = localMatrix(0, 0);
+}
+
+void BilinearTermBase::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const 
+{
+  std::vector<fullMatrix<double> > mv(npts);
+  get(ele,npts,GP,mv);
+  m.resize(mv[0].size1(), mv[0].size2());
+  m.setAll(0.);
+  double jac[3][3];
+  for (int k=0;k<npts;k++)
+  {
+    const double u = GP[k].pt[0]; const double v = GP[k].pt[1]; const double w = GP[k].pt[2];
+    const double weight = GP[k].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+    const double coeff=weight*detJ;
+    for (int i=0;i<mv[k].size1();++i)
+      for (int j=0;j<mv[k].size2();++j)
+        m(i,j)+=mv[k](i,j)*coeff;
+  }
+}
+
+IsotropicElasticTerm::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);
+  double C12 = FACT * nu / (1 - 2 * nu);
+  double C44 = (C11 - C12) / 2;
+  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_);
+}
+IsotropicElasticTerm::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);
+  double C12 = FACT * nu / (1 - 2 * nu);
+  double C44 = (C11 - C12) / 2;
+/*  FACT = E / (1 - nu * nu); // plane stress (plates)
+  C11  = FACT;
+  C12  = nu * FACT; 
+  C44 = (1. - nu) * .5 * FACT;*/
+  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;
+}
+
+void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+  if(ele->getParent()) ele = ele->getParent();
+  if(sym)
+  {
+    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.);
+    //std::cout << m.size1() << "  " << m.size2() << std::endl;
+    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.);
+    // Sum on Gauss Points i
+    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;// 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);
+      // gemm add the product to m so there is a sum on gauss' points here
+      m.gemm(BTH, B, weight * detJ, 1.);
+    }
+  }
+}
+
+void LagrangeMultiplierTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+  int nbFF1 = BilinearTerm<SVector3, double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+  int nbFF2 = BilinearTerm<SVector3, double>::space2.getNumKeys(ele); //nbVertices of boundary
+  double jac[3][3];
+  m.resize(nbFF1, nbFF2);
+  m.setAll(0.);
+  for(int i = 0; i < npts; i++)
+  {
+    double u = GP[i].pt[0]; double v = GP[i].pt[1]; 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>::ValType> Vals;
+    std::vector<TensorialTraits<double>::ValType> ValsT;
+    BilinearTerm<SVector3,double>::space1.f(ele, u, v, w, Vals);
+    BilinearTerm<SVector3,double>::space2.f(ele, u, v, w, ValsT);
+    for(int j = 0; j < nbFF1; j++)
+    {
+      for(int k = 0; k < nbFF2; k++)
+      {
+        m(j, k) += dot(Vals[j], _d) * ValsT[k] * weight * detJ;
+      }
+    }
+  }
+}
+
+void LagMultTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+  int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+  int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele); //nbVertices of boundary
+  double jac[3][3];
+  m.resize(nbFF1, nbFF2);
+  m.setAll(0.);
+  for(int i = 0; i < npts; i++) 
+  {
+    double u = GP[i].pt[0]; double v = GP[i].pt[1]; 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>::ValType> Vals;
+    std::vector<TensorialTraits<SVector3>::ValType> ValsT;
+    BilinearTerm<SVector3,SVector3>::space1.f(ele, u, v, w, Vals);
+    BilinearTerm<SVector3,SVector3>::space2.f(ele, u, v, w, ValsT);
+    for(int j = 0; j < nbFF1; j++)
+    {
+      for(int k = 0; k < nbFF2; k++)
+      {
+        m(j, k) += _eqfac * dot(Vals[j], ValsT[k]) * weight * detJ;
+      }
+    }
+  }
+}
+
+

Solver/terms.hpp

+//
+// C++ Template Implementations: terms
+//
+// Description:
+//
+//
+// Author:  <Eric Bechet>, (C) 2011
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "terms.h"
+
+template<class T2> void LinearTermBase<T2>::get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &vec) const
+{
+  std::vector<fullVector<T2> > vv;
+  vv.resize(npts);
+  get(ele,npts,GP,vv);
+  int nbFF=vv[0].size();
+  vec.resize(nbFF);
+  vec.setAll(T2());
+  double jac[3][3];
+  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);
+    for(int j = 0; j < nbFF; j++)
+    {
+      double contrib = weight * detJ;
+      vec(j) += contrib*vv[i](j);
+    }
+  }
+}
+
+template<class T2> void PlusTerm<T2>::get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &v) const 
+{
+  fullVector<T2> v2;
+  a->get(ele,npts,GP,v);
+  b->get(ele,npts,GP,v2);
+  for (int i=0;i<v2.size();++i) v(i)+=v2(i);
+}
+
+template<class T2> void ScalarTermConstant<T2>::get(MElement *ele, int npts, IntPt *GP, T2 &val) const
+{
+  double jac[3][3];
+  double eval = 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);
+    eval += weight * detJ;
+  }
+  val=cst * eval;
+}
+
+template<class T2> void ScalarTermConstant<T2>::get(MElement *ele, int npts, IntPt *GP, std::vector<T2> &vval) const
+{
+  for(int i = 0; i < npts; i++)
+  {
+    vval[i] = cst;
+  }
+}
+
+template<class T2> void ScalarTermConstant<T2>::get(MVertex *ver, T2 &val) const
+{
+    val = cst;
+}
+
+template<class T1> void LaplaceTerm<T1, T1>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+  int nbFF = BilinearTerm<T1, T1>::space1.getNumKeys(ele);
+  double jac[3][3];
+  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<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++)
+      {
+        double contrib = weight * detJ * dot(Grads[j], Grads[k]) * diffusivity;
+        m(j, k) += contrib;
+        if(j != k) m(k, j) += contrib;
+      }
+    }
+  }
+}
+
+template<class T1> void LoadTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const 
+{
+  if(ele->getParent()) ele = ele->getParent();
+  int nbFF = LinearTerm<T1>::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 TensorialTraits<T1>::ValType> Vals;
+    LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
+    SPoint3 p;
+    ele->pnt(u, v, w, p);
+    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;
+    }
+  }
+}
+
+template<class T2> BilinearTermContract<T2> operator |(const LinearTermBase<T2>& L1,const LinearTermBase<T2>& L2)
+{
+  return BilinearTermContract<T2>(L1,L2);
+}
+
+template<class T2> void BilinearTermContract<T2>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const 
+{
+  fullVector<T2> va;
+  fullVector<T2> vb;
+  a->get(ele,npts,GP,va);
+  b->get(ele,npts,GP,vb);
+  m.resize(va.size(), vb.size());
+  m.setAll(0.);
+  for (int i=0;i<va.size();++i)
+    for (int j=0;j<vb.size();++j)
+      m(i,j)=dot(va(i),vb(j));
+}
+
+
+template<class T2> void BilinearTermContractWithLaw<T2>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const 
+{
+  BilinearTermBase::get(ele,npts,GP,m);
+}
+  
+template<class T2> void BilinearTermContractWithLaw<T2>::get(MElement *ele, int npts, IntPt *GP, std::vector<fullMatrix<double> > &mv) const 
+{
+  std::vector<fullVector<T2> > va(npts);
+  std::vector<fullVector<T2> > vb(npts);
+  std::vector<typename TensorialTraits<T2>::TensProdType> tens(npts);
+  BilinearTermContract<T2>::a->get(ele,npts,GP,va);
+  BilinearTermContract<T2>::b->get(ele,npts,GP,vb);
+  c->get(ele,npts,GP,tens);
+  for (int k=0;k<npts;k++)
+  {
+    mv[k].resize(va[k].size(), vb[k].size());
+    for (int i=0;i<va[k].size();++i)
+      for (int j=0;j<vb[k].size();++j)
+        mv[k](i,j)=dot(va[k](i),tens[k]*vb[k](j));
+  }
+}
+
+template<class T2> PlusTerm<T2> LinearTermBase<T2>::operator +(const LinearTermBase<T2>& other)
+{
+  return PlusTerm<T2>(*this,other);
+}
+/*
+template<class T1> void GradTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<typename TensorialTraits<T1>::GradType > &vec) const 
+{
+  int nbFF = LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.getNumKeys(ele);
+  double jac[3][3];
+  vec.resize(nbFF);
+  vec.setAll(typename TensorialTraits<T1>::GradType());
+  for(int i = 0; i < npts; i++)
+  {
+    std::vector<typename TensorialTraits<T1>::GradType> Grads;
+    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);
+    LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.gradf(ele, u, v, w, Grads);
+    for(int j = 0; j < nbFF; j++)
+    {
+      double contrib = weight * detJ;
+      vec(j) += contrib*Grads[j];
+    }
+  }
+}
+*/
+template<class T1> void GradTerm<T1>::get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<typename TensorialTraits<T1>::GradType > > &vvec) const 
+{
+  int nbFF = LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.getNumKeys(ele);
+  for(int i = 0; i < npts; i++)
+  {
+    vvec[i].resize(nbFF);
+    std::vector<typename TensorialTraits<T1>::GradType> Grads;
+    const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
+    LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.gradf(ele, u, v, w, Grads);
+    for(int j = 0; j < nbFF; j++)
+    {
+      vvec[i](j) = Grads[j];
+    }
+  }
+}
+
+
+
+template<class T1> void LoadTermOnBorder<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const 
+{
+  MElement *elep;
+  if (ele->getParent()) elep = ele->getParent();
+  int nbFF = LinearTerm<T1>::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 TensorialTraits<T1>::ValType> Vals;
+    LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
+    SPoint3 p;
+    ele->pnt(u, v, w, p);
+    typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+    for(int j = 0; j < nbFF ; ++j){
+      m(j) += _eqfac * dot(Vals[j], load) * weight * detJ;
+    }
+  }
+}