no more mVertex* - related members in function spaces. Use Mpoint instead.

Solver/elasticitySolver.cpp

@@ -16,7 +16,7 @@
 #include "solverAlgorithms.h"
 #include "quadratureRules.h"
 #include "solverField.h"
-
+#include "MPoint.h"
 #if defined(HAVE_POST)
 #include "PView.h"
 #include "PViewData.h"
@@ -164,14 +164,11 @@ void elasticitySolver::solve()
   // give priority to fixations : when a dof is fixed, it cannot be
   // numbered afterwards
 
-
+  std::cout <<  "Dirichlet BC"<< std::endl;
   for (unsigned int i = 0; i < allDirichlet.size(); i++)
   {
     FilterDofComponent filter(allDirichlet[i]._comp);
-    if (allDirichlet[i].onWhat==BoundaryCondition::ON_VERTEX)
-      FixNodalDofs(*LagSpace,allDirichlet[i].g->vbegin(),allDirichlet[i].g->vend(),*pAssembler,allDirichlet[i]._f,filter);
-    else
-      FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
+    FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
   }
 
   // we number the dofs : when a dof is numbered, it cannot be numbered
@@ -185,14 +182,11 @@ void elasticitySolver::solve()
   // First build the force vector
 
   GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
-
+  std::cout <<  "Neumann BC"<< std::endl;
   for (unsigned int i = 0; i < allNeumann.size(); i++)
   {
     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
-      Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
+    Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
   }
 
 // bulk material law
@@ -279,7 +273,8 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
   for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
   {
     SVector3 val;
-    Field.f(*it,val);
+    MPoint p(*it);
+    Field.f(&p,0,0,0,val);
     std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
     data[(*it)->getNum()]=vec;
   }

Solver/function.h

@@ -51,6 +51,13 @@ public :
   inline bool somethingDependOnMe() {
     return !_dependOnMe.empty();
   }
+  inline int howManyDependOnMe() {
+    return _dependOnMe.size();
+  }
+  inline int howManyDoIDependOn() {
+    return _iDependOn.size();
+  }
+
 };
 
 // dataCache when the value is a double 

Solver/functionSpace.h

@@ -10,66 +10,40 @@
 #include "MElement.h"
 #include "dofManager.h"
 
-class SVoid{};
-
-class basisFunction{
- public:
-  virtual void f(MElement *ele, double u, double v, double w, double *s) = 0;
-  virtual void df(MElement *ele, double u, double v, double w, double *s) = 0;
-};
-
-class scalarPolynomialBasisFunction : public basisFunction{
- private:
-  polynomialBasis *_basis;
- public:
-  virtual void f(MElement *ele, double u, double v, double w, double *s);
-  virtual void df(MElement *ele, double u, double v, double w, double *s);
-};
-
-class vectorPolynomialBasisFunction : public basisFunction{
- private:
-  polynomialBasis *_basis[3];
- public:
-  virtual void f(MElement *ele, double u, double v, double w, double *s);
-  virtual void df(MElement *ele, double u, double v, double w, double *s);
-};
+//class SVoid{};
 
 template<class T> struct TensorialTraits
 {
   typedef T ValType;
   typedef T GradType[3];
-  typedef T HessType[3][3];
+/*  typedef T HessType[3][3];
   typedef SVoid DivType;
-  typedef SVoid CurlType;
+  typedef SVoid CurlType;*/
 };
 
 template<> struct TensorialTraits<double>
 {
   typedef double ValType;
   typedef SVector3 GradType;
-  typedef STensor3 HessType;
+/*  typedef STensor3 HessType;
   typedef SVoid DivType;
-  typedef SVoid CurlType;
+  typedef SVoid CurlType;*/
 };
 
 template<> struct TensorialTraits<SVector3>
 {
   typedef SVector3 ValType;
   typedef STensor3 GradType;
-  typedef SVoid HessType;
+/*  typedef SVoid HessType;
   typedef double DivType;
-  typedef SVector3 CurlType;
+  typedef SVector3 CurlType;*/
 };
 
-
-
 class FunctionSpaceBase
 {
  public:
   virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0;
-  virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
-  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
 };
 
 template<class T>
@@ -78,23 +52,10 @@ class FunctionSpace : public FunctionSpaceBase
  public:
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
-/*  typedef typename TensorialTraits<T>::HessType HessType;
-  typedef typename TensorialTraits<T>::DivType DivType;
-  typedef typename TensorialTraits<T>::CurlType CurlType;*/
- 
- public:
-  virtual void f(MVertex *ver, std::vector<ValType> &vals) {} // used for neumann BC
   virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0;
   virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)=0;
-//  virtual groupOfElements* getSupport()=0;// probablement inutile
-//  virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads, STensor3 &invjac)=0;// on passe le jacobien que l'on veut ...
-/* virtual int hessf(MElement *ele, double u, double v, double w,std::vector<HessType> &hesss);
-  virtual int divf(MElement *ele, double u, double v, double w,std::vector<DivType> &divs);
-  virtual int curlf(MElement *ele, double u, double v, double w,std::vector<CurlType> &curls);*/
   virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0;
-  virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
-  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
 };
 
 class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
@@ -102,14 +63,16 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
  public:
   typedef TensorialTraits<double>::ValType ValType;
   typedef TensorialTraits<double>::GradType GradType;
-  typedef TensorialTraits<double>::HessType HessType;
-  typedef TensorialTraits<double>::DivType DivType;
-  typedef TensorialTraits<double>::CurlType CurlType;
- protected:
-  std::vector<basisFunction*> basisFunctions;
 
+ protected:
   int _iField; // field number (used to build dof keys)
 
+ private:
+  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
+  {
+    keys.push_back(Dof(ver->getNum(), _iField));
+  }
+ 
  public:
   ScalarLagrangeFunctionSpace(int i=0):_iField(i) {}
   virtual int getId(void) const {return _iField;};
@@ -119,8 +82,8 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
     int curpos=vals.size();
     vals.resize(curpos+ndofs);
     ele->getShapeFunctions(u, v, w, &(vals[curpos]));
-  };
-  virtual void f(MVertex *ver, std::vector<ValType> &vals) {vals.push_back(1.0);} // used for neumann BC
+  }
+
   virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
   {
     int ndofs= ele->getNumVertices();
@@ -137,7 +100,8 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
       invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2],
       invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2]
                             ));
-  };
+  }
+
   virtual int getNumKeys(MElement *ele) {return ele->getNumVertices();}
 
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
@@ -147,11 +111,6 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
       for (int i=0;i<ndofs;++i)
         getKeys(ele->getVertex(i),keys);
   }
-  virtual int getNumKeys(MVertex *ver) { return 1;}
-  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
-  {
-    keys.push_back(Dof(ver->getNum(), _iField));
-  };
 };
 
 
@@ -160,9 +119,6 @@ template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T>
 public :
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
-  typedef typename TensorialTraits<T>::HessType HessType;
-  typedef typename TensorialTraits<T>::DivType DivType;
-  typedef typename TensorialTraits<T>::CurlType CurlType;
 protected :
   std::vector<T> multipliers;
   std::vector<int> comp;
@@ -187,19 +143,6 @@ public :
   }
 
   virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;}
-  virtual void f(MVertex *ver, std::vector<ValType> &vals)
-  {
-    std::vector<double> valsd;
-    ScalarFS->f(ver,valsd);
-    int nbdofs=valsd.size();
-    int nbcomp=comp.size();
-    int curpos=vals.size();
-    vals.reserve(curpos+nbcomp*nbdofs);
-    for (int j=0;j<nbcomp;++j)
-    {
-      for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]);
-    }
-  } // used for neumann BC
 
   virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
   {
@@ -257,28 +200,6 @@ public :
       }
     }
   }
-
-  virtual int getNumKeys(MVertex *ver) { return ScalarFS->getNumKeys(ver)*comp.size();}
-  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
-  {
-    int nk=ScalarFS->getNumKeys(ver);
-    std::vector<Dof> bufk;
-    bufk.reserve(nk);
-    ScalarFS->getKeys(ver,bufk);
-    int nbdofs=bufk.size();
-    int nbcomp=comp.size();
-    int curpos=keys.size();
-    keys.reserve(curpos+nbcomp*nbdofs);
-    for (int j=0;j<nbcomp;++j)
-    {
-      for (int i=0;i<nk;++i)
-      {
-        int i1,i2;
-        Dof::getTwoIntsFromType(bufk[i].getType(), i1,i2);
-        keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i1)));
-      }
-    }
-  };
 };
 
 class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
@@ -289,9 +210,6 @@ class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
   enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 };
   typedef TensorialTraits<SVector3>::ValType ValType;
   typedef TensorialTraits<SVector3>::GradType GradType;
-  typedef TensorialTraits<SVector3>::HessType HessType;
-  typedef TensorialTraits<SVector3>::DivType DivType;
-  typedef TensorialTraits<SVector3>::CurlType CurlType;
   VectorLagrangeFunctionSpace(int id) :
           ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
           SVector3(1.,0.,0.),VECTOR_X, SVector3(0.,1.,0.),VECTOR_Y,SVector3(0.,0.,1.),VECTOR_Z)
@@ -317,9 +235,6 @@ class CompositeFunctionSpace : public FunctionSpace<T>
  public: 
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
-  typedef typename TensorialTraits<T>::HessType HessType;
-  typedef typename TensorialTraits<T>::DivType DivType;
-  typedef typename TensorialTraits<T>::CurlType CurlType;
   typedef typename std::vector<FunctionSpace<T>* >::iterator iterFS;
  protected: 
  
@@ -347,13 +262,6 @@ class CompositeFunctionSpace : public FunctionSpace<T>
     for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
       delete (*it);
   }
-
-  virtual void f(MVertex *ver, std::vector<ValType> &vals)
-  {
-    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
-      (*it)->f(ver,vals);
-  }
-
   virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
   {
     for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
@@ -379,20 +287,6 @@ class CompositeFunctionSpace : public FunctionSpace<T>
     for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
       (*it)->getKeys(ele,keys);
   }
-  
-  virtual int getNumKeys(MVertex *ver)
-  {
-    int ndofs=0;
-    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
-      ndofs+=(*it)->getNumKeys(ver);
-    return ndofs;
-  }
-
-  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
-  {
-    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
-      (*it)->getKeys(ver,keys);
-  };
 };
 
 template<class T>
@@ -401,21 +295,15 @@ class xFemFunctionSpace : public FunctionSpace<T>
  public:
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
-  typedef typename TensorialTraits<T>::HessType HessType;
-  typedef typename TensorialTraits<T>::DivType DivType;
-  typedef typename TensorialTraits<T>::CurlType CurlType;
  private:
   FunctionSpace<T>* _spacebase;
 //  Function<double>* enrichment;
  public:
-  virtual void f(MVertex *ver, std::vector<ValType> &vals);
+
   virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
   virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
   virtual int getNumKeys(MElement *ele);
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
-  virtual int getNumKeys(MVertex *ver);
-  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
-
 };
 
 
@@ -424,25 +312,15 @@ class FilteredFunctionSpace : public FunctionSpace<T>
 {
   typedef typename TensorialTraits<T>::ValType ValType;
   typedef typename TensorialTraits<T>::GradType GradType;
-  typedef typename TensorialTraits<T>::HessType HessType;
-  typedef typename TensorialTraits<T>::DivType DivType;
-  typedef typename TensorialTraits<T>::CurlType CurlType;
  private:
   FunctionSpace<T>* _spacebase;
   F &_filter;
   
  public:
-  virtual void f(MVertex *ver, std::vector<ValType> &vals);
   virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
   virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
   virtual int getNumKeys(MElement *ele);
   virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
-  virtual int getNumKeys(MVertex *ver);
-  virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
 };
 
-
-
-
-
 #endif

Solver/solverAlgorithms.h

@@ -65,20 +65,6 @@ template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,Func
   }
 }
 
-template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
-{
-  fullVector<typename Assembler::dataMat> localVector;
-  std::vector<Dof> R;
-  for (Iterator it = itbegin;it!=itend; ++it)
-  {
-    MVertex *v = *it;
-    R.clear();
-    term.get(v,localVector);
-    space.getKeys(v,R);
-    assembler.assemble(R, localVector);
-  }
-}
-
 template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler)
 {
   fullVector<typename Assembler::dataMat> localVector;
@@ -177,20 +163,6 @@ template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e
   }
 }
 
-template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MVertex *v,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
-{
-  std::vector<Dof> R;
-  space.getKeys(v,R);
-  for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
-  {
-    Dof key=*itd;
-    if (filter(key))
-    {
-      assembler.fixDof(key, fct(v->x(),v->y(),v->z()));
-    }
-  }
-}
-
 template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
 {
   for (Iterator it=itbegin;it!=itend;++it)
@@ -211,5 +183,4 @@ template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &spac
   }
 }
 
-
 #endif// _SOLVERALGORITHMS_H_

Solver/solverField.h

@@ -21,7 +21,7 @@
 #include "functionSpace.h"
 
 template<class T>
-class SolverField : public FunctionSpace<T> // being able to use it instead of a real function space is interesting (nbkeys=1, dofs undefined (or could be defined by element-wise )
+class SolverField : public FunctionSpace<T> // being able to use it instead of a real function space is interesting (nbkeys=1, explicit keys/dofs undefined (or could be defined element-wise )
 {
  public:
   typedef typename TensorialTraits<T>::ValType ValType;
@@ -38,26 +38,6 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a
   virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) {Msg::Error("getKeys for SolverField should'nt be called");}
  public:
 
-  virtual void f(MVertex *ver, std::vector<ValType> &vals)
-  {
-    ValType val;
-    f(ver,val);
-    vals.push_back(val);
-  }
-
-  virtual void f(MVertex *ver, ValType &val)
-  {
-    std::vector<Dof> D;
-    std::vector<ValType> SFVals;
-    std::vector<double> DMVals;
-    fs->getKeys(ver,D);
-    dm->getDofValue(D,DMVals);
-    fs->f(ver,SFVals);
-    val=ValType();
-    for (int i=0;i<D.size();++i)
-      val+=SFVals[i]*DMVals[i];
-  }
-
   virtual void f(MElement *ele, double u, double v, double w, ValType &val)
   {
     std::vector<Dof> D;
@@ -99,9 +79,6 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a
   }
 };
 
-
-
-
 /*
 
 class Formulation

Solver/terms.h

@@ -44,7 +44,6 @@ 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 T1> class LinearTerm : public LinearTermBase
@@ -309,20 +308,6 @@ template<class T1> class LoadTerm : public LinearTerm<T1>
       }
     }
   }
-
-  virtual void get(MVertex *ver,fullVector<double> &m)
-  {
-    double nbFF=LinearTerm<T1>::space1.getNumKeys(ver);
-    double jac[3][3];
-    m.resize(nbFF);
-    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);
-    }
-  }
 };
 
 

contrib/arc/conge.dat

@@ -5,10 +5,9 @@ EdgeDisplacement 8 1 0
 EdgeDisplacement 8 2 0
 #NodeDisplacement 20 0 0
 #NodeDisplacement 20 1 0
-#NodeDisplacement 20 2 0
 #NodeDisplacement 21 1 0
-#NodeDisplacement 21 2 0
+#NodeForce 21 100e6 0 0
 EdgeForce 9 100e6 0 0
 #FaceForce 7  0 -1.0e10  0 
-#NodeForce 20 -1e7 0 0 
+#NodeForce 20 -1e8 0 0 
 

contrib/arc/mainElasticity.cpp

@@ -5,16 +5,106 @@
 #include "highlevel.h"
 #include "groupOfElements.h"
 #include <iterator>
+#include "function.h"
+#include "fullMatrix.h"
+/*
+class functionAdd : public function
+{
+ private:
+  std::vector<std::string> strvec;
+  class data : public dataCacheDouble
+  {
+   private:
+    std::vector<dataCacheDouble *> dcvec;
+   public:
+    data(const functionAdd * fm,dataCacheMap *m) :
+      dataCacheDouble(m->getNbEvaluationPoints(),1)
+    {
+      for (int i=0;i<fm->strvec.size();++i)
+        dcvec.push_back(&(m->get(fm->strvec[i],this)));
+    }
+    void _eval()
+    {
+      _value(0, 0)=0;
+      for (int i=0;i<dcvec.size();++i) { _value(0, 0) += (*(dcvec[i]))(0,0);}
+    }
+    ~data()
+    {
+    }
+  };
+ public:
+  dataCacheDouble *newDataCache(dataCacheMap *m)
+  {
+    return new data(this,m);
+  }
+  functionAdd(std::string _a,std::string _b) {strvec.push_back(_a);strvec.push_back(_b);}
+  void addNewTerm(std::string _a)  { strvec.push_back(_a);}
+};
+
+*/
+
 
 int main (int argc, char* argv[])
 {
-  
+
   if (argc != 2){
     printf("usage : elasticity input_file_name\n");
     return -1;
   }
+/*
+  fullMatrix<double> a(1,1);
+  a(0,0)=1.0;
+  fullMatrix<double> b(1,1);
+  b(0,0)=2.0;
+  fullMatrix<double> c(1,1);
+  c(0,0)=4.0;
+  fullMatrix<double> d(1,1);
+  d(0,0)=8.0;
+  fullMatrix<double> e(1,1);
+  e(0,0)=16.0;
+
+  std::cout << argc << std::endl;
+  functionConstant fc_a(&a);
+  functionConstant fc_b(&b);
+  functionConstant fc_c(&c);
+  functionConstant fc_d(&d);
+  functionConstant fc_e(&e);
+
+  functionAdd fa_ab(fc_a.getName(),fc_b.getName());
+  functionAdd fa_cd(fc_c.getName(),fc_d.getName());
+  functionAdd fa_abcd("a+b","c+d");
+
+  function::add("a+b", &fa_ab );
+  function::add("c+d", &fa_cd);
+  function::add("a+b+c+d", &fa_abcd);
+
+  dataCacheMap m(1);
+  std::cout << "start" << std::endl;
+  dataCacheDouble &dc_a=m.get(fc_a.getName());
+  dataCacheDouble &dc_abcd=m.get("a+b+c+d");
+  fa_abcd.addNewTerm(fc_e.getName());
+  std::cout << "nb depsabcd=" << dc_abcd.howManyDoIDependOn() << std::endl;
+  std::cout << "nb depsa=" << dc_a.howManyDependOnMe() << std::endl;
+  std::cout << "a" << std::endl;
+  std::cout << dc_a(0,0) << std::endl;
+  std::cout << "a+b+c+d" << std::endl;
+  std::cout << dc_abcd(0,0) << std::endl;
+  std::cout << "dca.set(b)" << std::endl;
+  dc_a.set(b);
+  std::cout << "a+b+c+d" << std::endl;
+  std::cout << dc_abcd(0,0) << std::endl;
+*/
+
+/*
+  functionMult fm("axbxcxd","axbxcxd");
+  dataCacheDouble *res;
+  res=fm.newDataCache(&m);
+  std::cout << "*res" << std::endl;
+  std::cout << (*res)(0,0) << std::endl;
+*/
+
+//  return(0);
 
-  
   GmshInitialize(argc, argv);
   // globals are still present in Gmsh
 
@@ -38,7 +128,7 @@ int main (int argc, char* argv[])
 
   // stop gmsh
   GmshFinalize();
-  
+
 }
 
 
@@ -60,13 +150,13 @@ int main (int argc, char* argv[])
   ScalarLagrangeFunctionSpace L(100);
   std::cout << L.getNumKeys(e) << "fonctions de formes L" << std::endl;
   L.getKeys(e,dofs);
-  for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl; 
+  for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl;
   dofs.clear();
   L.f(e,0.1,0.1,0,vals);
   L.gradf(e,0.1,0.1,0,grads);
   std::copy(vals.begin(),vals.end(),output); std::cout << std::endl;
   for (std::vector<SVector3>::iterator it=grads.begin();it!=grads.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; }
-  
+
   VectorLagrangeFunctionSpace L1(100,VectorLagrangeFunctionSpace::VECTOR_X);
   VectorLagrangeFunctionSpace L2(100,VectorLagrangeFunctionSpace::VECTOR_Y);
   std::cout << L2.getNumKeys(e) << "fonctions de formes L2" << std::endl;
@@ -81,7 +171,7 @@ int main (int argc, char* argv[])
   std::cout << P123.getNumKeys(e) << "fonctions de formes P123" << std::endl;
   P123.getKeys(e,dofs);
   std::cout << dofs.size() << std::endl;
-  for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl; 
+  for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl;
 
   vals2.clear();
   grads2.clear();
@@ -99,7 +189,7 @@ int main (int argc, char* argv[])
 
   FormBilinear<TermBilinearMecaNL,ScalarLagrangeFunctionSpace,ScalarLagrangeFunctionSpace > fnl(L,L);
   fnl.func();
-  
+
 */