Messed up with function spaces / working example in linear elasiticity

Solver/CMakeLists.txt

@@ -21,6 +21,7 @@ set(SRC
   function.cpp
   functionDerivator.cpp
   luaFunction.cpp
+  functionSpace.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 

Solver/SElement.cpp

@@ -42,7 +42,7 @@ void SElement::gradNodalFunctions (double u, double v, double w, double invjac[3
     _e->getShapeFunctions(u, v, w, sf);
     // FIXME : re-use sf for computing coordinates
     _e->pnt(u, v, w, p);
-    double E = (*_enrichement_s)(p.x(), p.y(), p.z());
+    double E = (*_enrichement)(p.x(), p.y(), p.z());
     double dEdx, dEdy, dEdz;
     _enrichement_s->gradient(p.x(), p.y(), p.z(), dEdx, dEdy, dEdz);
     for (int i = 0; i < N; i++){
@@ -62,7 +62,7 @@ void SElement::nodalFunctions (double u, double v, double w, double s[],
     SPoint3 p;
     // FIXME : re-use s for computing coordinates
     _e->pnt(u, v, w, p);
-    double E = (*_enrichement_s)(p.x(), p.y(), p.z());
+    double E = (*_enrichement)(p.x(), p.y(), p.z());
     for (int i = 0; i < N; i++){
       s[i] *= E;
     }

Solver/dofManager.h

@@ -23,11 +23,11 @@ class Dof{
   Dof(long int entity, int type) : _entity(entity), _type(type) {}
   inline long int getEntity() const { return _entity; }
   inline int getType() const { return _type; }
-  static int createTypeWithTwoInts(int i1, int i2)
+  inline static int createTypeWithTwoInts(int i1, int i2)
   {
     return i1 + 10000 * i2;
   }
-  static void getTwoIntsFromType(int t, int &i1, int &i2)
+  inline static void getTwoIntsFromType(int t, int &i1, int &i2)
   {
     i1 = t % 10000;
     i2 = t / 10000;
@@ -50,18 +50,18 @@ template<class T> struct dofTraits
   typedef T MatType;
 };
 
-template<> struct dofTraits<fullVector<double> >
+template<class T> struct dofTraits<fullVector<T> >
 {
-  typedef fullVector<double> VecType;
-  typedef fullMatrix<double> MatType;
+  typedef fullVector<T> VecType;
+  typedef fullMatrix<T> MatType;
 };
-
+/*
 template<> struct dofTraits<fullVector<std::complex<double> > >
 {
   typedef fullVector<std::complex<double> > VecType;
   typedef fullMatrix<std::complex<double> > MatType;
 };
-
+*/
 template<class T>
 class DofAffineConstraint{
  public:

Solver/elasticitySolver.cpp

@@ -13,6 +13,8 @@
 #include "linearSystemPETSc.h"
 #include "linearSystemGMM.h"
 #include "Numeric.h"
+#include "functionSpace.h"
+#include "terms.h"
 
 #if defined(HAVE_POST)
 #include "PView.h"
@@ -321,3 +323,171 @@ void elasticitySolver::solve()
   // solving
   lsys->systemSolve();
 }
+
+
+
+
+void MyelasticitySolver::solve()
+{
+#if defined(HAVE_TAUCS)
+  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+  linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#else
+  linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+  lsys->setNoisy(2);
+#endif
+  pAssembler = new dofManager<double>(lsys);
+
+  std::map<int, std::vector<GEntity*> > groups[4];
+  pModel->getPhysicalGroups(groups);
+
+  // we first do all fixations. the behavior of the dofManager is to 
+  // give priority to fixations : when a dof is fixed, it cannot be
+  // numbered afterwards
+
+  for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
+       it != nodalDisplacements.end(); ++it){
+    MVertex *v = pModel->getMeshVertexByTag(it->first.first);
+    pAssembler->fixVertex(v, it->first.second, _tag, it->second);
+    printf("-- Fixing node %3d comp %1d to %8.5f\n",
+           it->first.first, it->first.second, it->second);
+  }
+
+  for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
+       it != edgeDisplacements.end(); ++it){
+    DummyfemTerm El(pModel);
+    El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag,
+                        simpleFunction<double>(it->second), *pAssembler);
+    printf("-- Fixing edge %3d comp %1d to %8.5f\n",
+           it->first.first, it->first.second, it->second);
+  }
+
+  for (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
+       it != faceDisplacements.end(); ++it){
+    DummyfemTerm El(pModel);
+    El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag,
+                        simpleFunction<double>(it->second), *pAssembler);
+    printf("-- Fixing face %3d comp %1d to %8.5f\n",
+           it->first.first, it->first.second, it->second);
+  }
+
+  // we number the nodes : when a node is numbered, it cannot be numbered
+  // again with another number. so we loop over elements
+  for (unsigned int i = 0; i < elasticFields.size(); ++i){
+    groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+    for ( ; it != elasticFields[i].g->end(); ++it){
+      SElement se(*it);
+      for (int j = 0; j < se.getNumNodalShapeFunctions(); ++j){
+        pAssembler->numberVertex(se.getVertex(j), 0, elasticFields[i]._tag);
+        pAssembler->numberVertex(se.getVertex(j), 1, elasticFields[i]._tag);
+        pAssembler->numberVertex(se.getVertex(j), 2, elasticFields[i]._tag);	
+      }
+    }
+  }
+
+  // Now we start the assembly process
+  // First build the force vector
+  // Nodes at geometrical vertices
+  for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
+       it != nodalForces.end(); ++it){
+    int iVertex = it->first;
+    SVector3 f = it->second;
+    std::vector<GEntity*> ent = groups[0][iVertex];
+    for (unsigned int i = 0; i < ent.size(); i++){
+      pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
+      pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
+      pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
+      printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
+             ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
+    }
+  }
+
+  // line forces
+  for (std::map<int, SVector3 >::iterator it = lineForces.begin();
+       it != lineForces.end(); ++it){
+    DummyfemTerm El(pModel);
+    int iEdge = it->first;
+    SVector3 f = it->second;
+    El.neumannNodalBC(iEdge, 1, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
+    El.neumannNodalBC(iEdge, 1, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
+    El.neumannNodalBC(iEdge, 1, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+    printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
+  }
+
+  // face forces
+  for (std::map<int, SVector3 >::iterator it = faceForces.begin();
+       it != faceForces.end(); ++it){
+    DummyfemTerm El(pModel);
+    int iFace = it->first;
+    SVector3 f = it->second;
+    El.neumannNodalBC(iFace, 2, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
+    El.neumannNodalBC(iFace, 2, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
+    El.neumannNodalBC(iFace, 2, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+    printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
+  }
+
+  // volume forces
+  for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
+       it != volumeForces.end(); ++it){
+    DummyfemTerm El(pModel);
+    int iVolume = it->first;
+    SVector3 f = it->second;
+//    El.setVector(f);
+    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]);
+    }
+  }
+
+  // assembling the stifness matrix
+  for (unsigned int i = 0; i < elasticFields.size(); i++){
+    SElement::setShapeEnrichement(elasticFields[i]._enrichment);
+    for (unsigned int j = 0; j < elasticFields.size(); j++){
+      elasticityTerm El(0, elasticFields[i]._E, elasticFields[i]._nu,
+                        elasticFields[i]._tag, elasticFields[j]._tag);
+      SElement::setTestEnrichement(elasticFields[j]._enrichment);
+      //      printf("%d %d\n",elasticFields[i]._tag,elasticFields[j]._tag);
+      El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);      
+    }
+  }
+
+/*  VectorLagrangeFunctionSpace L1(Dof::createTypeWithTwoInts(0,_tag),VectorLagrangeFunctionSpace::VECTOR_X);
+  VectorLagrangeFunctionSpace L2(Dof::createTypeWithTwoInts(1,_tag),VectorLagrangeFunctionSpace::VECTOR_Y);
+  VectorLagrangeFunctionSpace L3(Dof::createTypeWithTwoInts(2,_tag),VectorLagrangeFunctionSpace::VECTOR_Z);
+  CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType> P123(L1,L2,L3);*/
+  VectorLagrangeFunctionSpace P123(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+
+  for (unsigned int i = 0; i < elasticFields.size(); i++)
+  {
+    DummyfemTerm El(pModel);
+//    ElasticTerm<CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType>,CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType> > Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
+    ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
+    fullMatrix<double> localMatrix(12,12);
+    std::vector<Dof> R;R.reserve(100);
+    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);
+      El.addToMatrix(*pAssembler,localMatrix,R);
+    }
+  }
+
+  printf("-- done assembling!\n");
+  //  for (int i=0;i<pAssembler->sizeOfR();i++){
+    //    printf("%g ",lsys->getFromRightHandSide(i));
+  //  }
+  //  printf("\n");
+
+  // solving
+  lsys->systemSolve();
+}
+
+

Solver/elasticitySolver.h

@@ -63,7 +63,7 @@ class elasticitySolver{
   void setMesh(const std::string &meshFileName);
   virtual void solve();
   void readInputFile(const std::string &meshFileName);
-  PView *buildDisplacementView(const std::string &postFileName);
+  virtual PView *buildDisplacementView(const std::string &postFileName);
   // PView *buildVonMisesView(const std::string &postFileName);
   // std::pair<PView *, PView*> buildErrorEstimateView
   //   (const std::string &errorFileName, double, int);
@@ -71,4 +71,16 @@ class elasticitySolver{
   //   (const std::string &errorFileName, const elasticityData &ref, double, int);
 };
 
+
+
+// another elastic solver ...
+class MyelasticitySolver : public elasticitySolver
+{
+ public:
+  MyelasticitySolver(int tag) : elasticitySolver(tag) {}
+  virtual void solve();
+};
+
+
+
 #endif

Solver/femTerm.h

@@ -179,4 +179,42 @@ class femTerm {
   }
 };
 
+
+
+class DummyfemTerm : public femTerm<double>
+{
+ public:
+  typedef dofTraits<double>::VecType dataVec;
+  typedef dofTraits<double>::MatType dataMat;
+  DummyfemTerm(GModel *gm) : femTerm<double>(gm) {}
+  virtual ~DummyfemTerm() {}
+ private : // i dont want to mess with this anymore
+  virtual int sizeOfC(SElement *se) const {return 0;}
+  virtual int sizeOfR(SElement *se) const {return 0;}
+  virtual Dof getLocalDofR(SElement *se, int iRow) const {return Dof(0,0);}
+  virtual Dof getLocalDofC(SElement *se, int iCol) const {return Dof(0,0);}
+  virtual void elementMatrix(SElement *se, fullMatrix<dataMat> &m) const {m.scale(0.);}
+  virtual void elementVector(SElement *se, fullVector<dataVec> &m) const {m.scale(0.);}
+ public:
+  void addToMatrix(dofManager<dataVec> &dm,
+                   fullMatrix<dataMat> &localMatrix,
+                   std::vector<Dof> R,std::vector<Dof> C) const
+  {
+      dm.assemble(R, C, localMatrix);
+  }
+
+
+  void addToMatrix(dofManager<dataVec> &dm,
+                   fullMatrix<dataMat> &localMatrix,
+                   std::vector<Dof> R) const // symmetric version.
+  {
+      dm.assemble(R, localMatrix);
+  }
+
+};
+
+
+
+
+
 #endif

Solver/functionSpace.cpp

+//
+// C++ Implementation: functionSpace
+//
+// Description: 
+//
+//
+// Author:  <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include "functionSpace.h"
+
+const SVector3 VectorLagrangeFunctionSpace::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};

Solver/functionSpace.h

@@ -2,12 +2,14 @@
 #define _FUNCTION_SPACE_H_
 
 #include "SVector3.h"
+#include "STensor3.h"
 #include <vector>
 #include <iterator>
 #include "Numeric.h"
+#include "MElement.h"
+#include "dofManager.h"
 
-
-class STensor3{};
+//class STensor3{};
 class SVoid{};
 
 class basisFunction{
@@ -71,6 +73,7 @@ class FunctionSpace {
  public:
   virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0;
   virtual int 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);
@@ -112,7 +115,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
   virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
   {
     int ndofs= ele->getNumVertices();
-    grads.reserve(grads.size()+ndofs);
+//    grads.reserve(grads.size()+ndofs);
     double gradsuvw[256][3];
     ele->getGradShapeFunctions(u, v, w, gradsuvw);
     double jac[3][3];
@@ -136,12 +139,13 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
   virtual int getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
   {
       int ndofs= ele->getNumVertices();
-      keys.reserve(keys.size()+ndofs);
+//      keys.reserve(keys.size()+ndofs);
       for (int i=0;i<ndofs;++i)
         keys.push_back(getLocalDof(ele,i));
   }
 };
 
+/*
 template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T> // scalarFS * const vector (avec vecteur non const, peut etre utilise pour xfem directement
 {
 public :
@@ -163,41 +167,160 @@ public :
     int nbdofs=valsd.size();
     for (int i=0;i<nbdofs;++i) vals.push_back(multiplier*valsd[i]);
   }
-  virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads){};
+  
+  virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+  {
+    std::vector<SVector3> gradsd;
+    ScalarFS->gradf(ele,u,v,w,gradsd);
+    int nbdofs=gradsd.size();
+    GradType val;
+    for (int i=0;i<nbdofs;++i)
+    {
+      tensprod(multiplier,gradsd[i],val);
+      grads.push_back(val);
+    }
+  };
   virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele);}
   virtual int getKeys(MElement *ele, Dof *keys){ ScalarFS->getKeys(ele,keys);}
   virtual int getKeys(MElement *ele, std::vector<Dof> &keys) {ScalarFS->getKeys(ele,keys);}
 };
+*/
 
+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;
+  FunctionSpace<double> *ScalarFS;
+public : 
+  template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult, int comp_): ScalarFS(new T2(SFS))
+  {
+    multipliers.push_back(mult);comp.push_back(comp_);
+  }
+
+  template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult1, int comp1_, 
+                          const T& mult2, int comp2_): ScalarFS(new T2(SFS)) 
+  {
+    multipliers.push_back(mult1);multipliers.push_back(mult2);comp.push_back(comp1_);comp.push_back(comp2_);
+  }
+
+  template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult1, int comp1_, 
+                          const T& mult2, int comp2_,const T& mult3, int comp3_): ScalarFS(new T2(SFS)) 
+  {
+    multipliers.push_back(mult1);multipliers.push_back(mult2);multipliers.push_back(mult3);
+    comp.push_back(comp1_);comp.push_back(comp2_);comp.push_back(comp3_);
+  }
 
-class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3> // it is a scalar lagrange times a constant vector.
+  virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;}
+  virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+  {
+    std::vector<double> valsd;
+    ScalarFS->f(ele,u,v,w,valsd);
+    int nbdofs=valsd.size();
+    int nbcomp=comp.size();
+    for (int j=0;j<nbcomp;++j)
+    {
+      for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]);
+    }
+  }
+  
+  virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+  {
+    std::vector<SVector3> gradsd;
+    ScalarFS->gradf(ele,u,v,w,gradsd);
+    int nbdofs=gradsd.size();
+    int nbcomp=comp.size();
+    GradType val;
+    for (int j=0;j<nbcomp;++j)
+    {
+      for (int i=0;i<nbdofs;++i)
+      {
+        tensprod(multipliers[j],gradsd[i],val);
+        grads.push_back(val);
+      }
+    }
+  };
+  virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele)*comp.size();}
+  virtual int getKeys(MElement *ele, Dof *keys)
+  { 
+    int nbcomp=comp.size();
+    int nk=ScalarFS->getNumKeys(ele);
+    Dof *kptr=keys;
+    ScalarFS->getKeys(ele,kptr);
+    kptr+=nk;
+    for (int j=1;j<nbcomp;++j)
+    {
+      for (int i=0;i<nk;++i)
+        kptr[i]=kptr[i-nk];
+      kptr+=nk;
+    }
+    kptr=keys;
+    for (int j=0;j<nbcomp;++j)
+    {
+      for (int i=0;i<nk;++i)
+      {
+        int i1,i2;
+        Dof::getTwoIntsFromType(kptr[i].getType(), i1,i2);
+        kptr[i]=Dof(kptr[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i2));
+      }
+      kptr+=nk;
+    }
+  }
+  virtual int getKeys(MElement *ele, std::vector<Dof> &keys)
+  {
+    int nk=ScalarFS->getNumKeys(ele);
+    std::vector<Dof> bufk;
+    bufk.reserve(nk);
+    ScalarFS->getKeys(ele,bufk);
+    int nbcomp=comp.size();
+    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],i2)));
+      }
+    }
+  }
+};
+
+class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
 {
  public:
   enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 };
+  static const SVector3 BasisVectors[3];
+
   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;
- protected:
-//  Along direction;
-  public:
-  VectorLagrangeFunctionSpace(int id,Along t) : ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),SVector3(0.,0.,0.) )//,direction(t)
-  { 
-    multiplier[t]=1.;
-  }
-//   virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals);
-//   virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
-//   virtual int getNumKeys(MElement *ele);
-//   virtual int getKeys(MElement *ele, Dof *keys);
-//   virtual int getKeys(MElement *ele, std::vector<Dof> &keys);
+  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)
+  {}
+  VectorLagrangeFunctionSpace(int id,Along comp1) :
+          ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
+          BasisVectors[comp1],comp1)
+  {}
+  VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2) :
+          ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
+          BasisVectors[comp1],comp1, BasisVectors[comp2],comp2)
+  {}
+  VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2, Along comp3) :
+          ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
+          BasisVectors[comp1],comp3, BasisVectors[comp2],comp2, BasisVectors[comp3],comp3)
+  {}
 };
 
 
-
-
-
-
 template<class T>
 class CompositeFunctionSpace : public FunctionSpace<T>
 {
@@ -235,8 +358,18 @@ class CompositeFunctionSpace : public FunctionSpace<T>
       delete (*it);
   }
 
-  virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals) {}
-  virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) {}
+  virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
+  {
+    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+      (*it)->f(ele,u,v,w,vals);
+  }
+
+  virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+  {
+    for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+      (*it)->gradf(ele,u,v,w,grads);
+  }
+
   virtual int getNumKeys(MElement *ele)
   {
     int ndofs=0;
@@ -244,6 +377,7 @@ class CompositeFunctionSpace : public FunctionSpace<T>
       ndofs+=(*it)->getNumKeys(ele);
     return ndofs;
   }
+
   virtual int getKeys(MElement *ele, Dof *keys)
   {
     Dof *kptr=keys;
@@ -280,4 +414,28 @@ class xFemFunctionSpace : public FunctionSpace<T>
 };
 
 
+template<class T,class F>
+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:
+  ValType f(MElement *ele, double u, double v, double w);
+  GradType gradf(MElement *ele, double u, double v, double w);
+  int getNumKeys(MElement *ele);
+  void getKeys(MElement *ele, Dof *keys);
+  void getKeys(MElement *ele, std::vector<Dof> &keys);
+};
+
+
+
+
+
 #endif