more generic algorithms

Solver/dofManager.h

@@ -74,10 +74,10 @@ class DofAffineConstraint{
 // fullVecor<double>, ...
 template <class T>
 class dofManager{
- private:
+ public:
   typedef typename dofTraits<T>::VecType dataVec;
   typedef typename dofTraits<T>::MatType dataMat;
-
+ private:
   // general affine constraint on sub-blocks, treated by adding
   // equations:
   //   dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
@@ -102,6 +102,10 @@ class dofManager{
    
  public:
   dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
+  inline void fixDof(Dof key, const dataVec &value)
+  {
+    fixed[key] = value;
+  }
   inline void fixDof(long int ent, int type, const dataVec &value)
   {
     fixed[Dof(ent, type)] = value;
@@ -110,9 +114,8 @@ class dofManager{
   {
     fixDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
   }
-  inline void numberDof(long int ent, int type)
+  inline void numberDof(Dof key)
   {
-    Dof key(ent, type);
     if(fixed.find(key) != fixed.end()) return;
     // if (constraints.find(key) != constraints.end()) return;
     std::map<Dof, int> :: iterator it = unknown.find(key);
@@ -121,6 +124,10 @@ class dofManager{
       unknown[key] = size;
     }
   }
+  inline void numberDof(long int ent, int type)
+  {
+    numberDof(Dof(ent, type));
+  }
   inline void numberVertex(MVertex*v, int iComp, int iField)
   {
     numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));

Solver/elasticitySolver.cpp

@@ -16,6 +16,7 @@
 #include "functionSpace.h"
 #include "terms.h"
 #include "solverAlgorithms.h"
+#include "quadratureRules.h"
 
 #if defined(HAVE_POST)
 #include "PView.h"
@@ -143,12 +144,24 @@ void elasticitySolver::readInputFile(const std::string &fn)
       int edge, comp;
       if(fscanf(f, "%d %d %lf", &edge, &comp, &val) != 3) return;
       edgeDisplacements[ std::make_pair(edge, comp) ] = val;
+      dirichletBC diri;
+      diri.g = new groupOfElements (1, edge);
+      diri._f= simpleFunction<double>(val);
+      diri._comp=comp;
+      diri._tag=edge;
+      edgeDiri.push_back(diri);
     }
     else if (!strcmp(what, "FaceDisplacement")){
       double val;
       int face, comp;
       if(fscanf(f, "%d %d %lf", &face, &comp, &val) != 3) return;
       faceDisplacements[ std::make_pair(face, comp) ] = val;
+      dirichletBC diri;
+      diri.g = new groupOfElements (2, face);
+      diri._f= simpleFunction<double>(val);
+      diri._comp=comp;
+      diri._tag=face;
+      faceDiri.push_back(diri);    
     }
     else if (!strcmp(what, "NodalForce")){
       double val1, val2, val3;
@@ -158,22 +171,37 @@ void elasticitySolver::readInputFile(const std::string &fn)
     }
     else if (!strcmp(what, "LineForce")){
       double val1, val2, val3;
-      int node;
-      if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
+      int edge;
+      if(fscanf(f, "%d %lf %lf %lf", &edge, &val1, &val2, &val3) != 4) return;
       //printf("%d %lf %lf %lf\n", node, val1, val2, val3);
-      lineForces[node] = SVector3(val1, val2, val3);    
+      lineForces[edge] = SVector3(val1, val2, val3); 
+      neumannBC neu;
+      neu.g = new groupOfElements (1, edge);
+      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._tag=edge;
+      edgeNeu.push_back(neu);
     }
     else if (!strcmp(what, "FaceForce")){
       double val1, val2, val3;
       int face;
       if(fscanf(f, "%d %lf %lf %lf", &face, &val1, &val2, &val3) != 4) return;
       faceForces[face] = SVector3(val1, val2, val3);
+      neumannBC neu;
+      neu.g = new groupOfElements (2, face);
+      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._tag=face;
+      edgeNeu.push_back(neu);
     }
     else if (!strcmp(what, "VolumeForce")){
       double val1, val2, val3;
       int volume;
       if(fscanf(f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3) != 4) return;
       volumeForces[volume] = SVector3(val1, val2, val3);
+      neumannBC neu;
+      neu.g = new groupOfElements (3, volume);
+      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._tag=volume;
+      edgeNeu.push_back(neu);  
     }
     else if (!strcmp(what, "MeshFile")){
       char name[245];
@@ -349,6 +377,7 @@ void MyelasticitySolver::solve()
   // 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
+  VectorLagrangeFunctionSpace P123(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
 
   for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
        it != nodalDisplacements.end(); ++it){
@@ -358,36 +387,24 @@ void MyelasticitySolver::solve()
            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 (unsigned int i = 0; i < edgeDiri.size(); i++)
+  {
+    FixNodalDofs(P123,edgeDiri[i].g->begin(),edgeDiri[i].g->end(),*pAssembler,edgeDiri[i]._f,FilterDofComponent(edgeDiri[i]._comp));
+    printf("-- Fixing edge %3d comp %3d \n", edgeDiri[i]._tag, edgeDiri[i]._comp);
   }
 
-  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);
+  for (unsigned int i = 0; i < faceDiri.size(); i++)
+  {
+    FixNodalDofs(P123,faceDiri[i].g->begin(),faceDiri[i].g->end(),*pAssembler,faceDiri[i]._f,FilterDofComponent(faceDiri[i]._comp));
+    printf("-- Fixing face %3d comp %3d \n", faceDiri[i]._tag, faceDiri[i]._comp);
   }
 
   // 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);	
-      }
-    }
+  for (unsigned int i = 0; i < elasticFields.size(); ++i)
+  {
+    NumberDofs(P123, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
   }
 
   // Now we start the assembly process
@@ -407,91 +424,37 @@ void MyelasticitySolver::solve()
     }
   }
 
-  VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
 
-  //line forces
-  for (std::map<int, SVector3 >::iterator it =lineForces.begin();
-       it != lineForces.end(); ++it)
-  {
-    int iEdge = it->first;
-    SVector3 f = it->second;
-    LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
-    printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
-    int dim=1;
-    std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iEdge);
-    if (itg == groups[dim].end())  {printf(" Nonexistent edge\n");break;}
-    for (unsigned int i = 0; i < itg->second.size(); ++i)
-    {
-      GEntity *ge = itg->second[i];
-      for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
+  GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
+
+  for (unsigned int i = 0; i < edgeNeu.size(); i++)
   {
-        MElement *e = ge->getMeshElement(j);
-        Assemble(Lterm,P123,e,*pAssembler);
-      }
+    LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,edgeNeu[i]._f);
+    Assemble(Lterm,P123,edgeNeu[i].g->begin(),edgeNeu[i].g->end(),Integ_Boundary,*pAssembler);
+    printf("-- Force on edge %3d \n", edgeNeu[i]._tag);
   }
-  }
-
 
-//face forces
-  for (std::map<int, SVector3 >::iterator it = faceForces.begin();
-       it != faceForces.end(); ++it)
+  for (unsigned int i = 0; i < faceNeu.size(); i++)
   {
-    int iFace = it->first;
-    SVector3 f = it->second;
-    LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
-    printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
-    int dim=2;
-    std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iFace);
-    if (itg == groups[dim].end())  {printf(" Nonexistent face\n");break;}
-    for (unsigned int i = 0; i < itg->second.size(); ++i)
-    {
-      GEntity *ge = itg->second[i];
-      for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
-      {
-        MElement *e = ge->getMeshElement(j);
-        Assemble(Lterm,P123,e,*pAssembler);
-      }
+    LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,faceNeu[i]._f);
+    Assemble(Lterm,P123,faceNeu[i].g->begin(),faceNeu[i].g->end(),Integ_Boundary,*pAssembler);
+    printf("-- Force on face %3d \n", faceNeu[i]._tag);
   }
-  }
-
 
-//volume forces
-  for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
-       it != volumeForces.end(); ++it)
+  for (unsigned int i = 0; i < volumeNeu.size(); i++)
   {
-    int iVolume = it->first;
-    SVector3 f = it->second;
-    LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
-    printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
-    int dim=3;
-    std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iVolume);
-    if (itg == groups[dim].end()) {printf(" Nonexistent volume\n");break;}
-    for (unsigned int i = 0; i < itg->second.size(); ++i)
-    {
-      GEntity *ge = itg->second[i];
-      for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
-      {
-        MElement *e = ge->getMeshElement(j);
-        Assemble(Lterm,P123,e,*pAssembler);
-      }
+    LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,volumeNeu[i]._f);
+    Assemble(Lterm,P123,volumeNeu[i].g->begin(),volumeNeu[i].g->end(),Integ_Boundary,*pAssembler);
+    printf("-- Force on volume %3d \n", volumeNeu[i]._tag);
   }
-  }
-
 
+  GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
   for (unsigned int i = 0; i < elasticFields.size(); i++)
   {
     ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
-    Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),*pAssembler);
+    Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
   }
 
-/*
-    for (int i=0;i<pAssembler->sizeOfR();i++)
-    {
-        if (fabs(lsys->getFromRightHandSide(i))>0.000001)
-        printf("%g ",lsys->getFromRightHandSide(i));
-    }
-    printf("\n");
-*/
   printf("-- done assembling!\n");
   lsys->systemSolve();
   printf("-- done solving!\n");

Solver/elasticitySolver.h

@@ -21,7 +21,22 @@ struct elasticField {
   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){}
+  elasticField () : g(0), _enrichment(0),_tag(0){}
+};
+
+struct dirichletBC {
+  int _tag; // tag for the dofManager
+  int _comp; // component
+  groupOfElements *g; // support for this BC
+  simpleFunction<double> _f;
+  dirichletBC () : g(0),_comp(0),_tag(0){}
+};
+
+struct neumannBC {
+  int _tag; // tag for the dofManager
+  groupOfElements *g; // support for this BC
+  simpleFunction<SVector3> _f;
+  neumannBC () : g(0),_tag(0),_f(SVector3(0,0,0)){}
 };
 
 // an elastic solver ...
@@ -36,16 +51,22 @@ class elasticitySolver{
   std::map<int, SVector3> nodalForces;
   // imposed line forces
   std::map<int, SVector3> lineForces;
+  std::vector<neumannBC> edgeNeu;
   // imposed face forces
   std::map<int, SVector3> faceForces;
+  std::vector<neumannBC> faceNeu;
   // imposed volume forces
   std::map<int, SVector3> volumeForces;
+  std::vector<neumannBC> volumeNeu;
   // imposed nodal displacements
   std::map<std::pair<int,int>, double> nodalDisplacements;
   // imposed edge displacements
   std::map<std::pair<int,int>, double> edgeDisplacements;
+  std::vector<dirichletBC> edgeDiri;
   // imposed face displacements
   std::map<std::pair<int,int>, double> faceDisplacements;
+  std::vector<dirichletBC> faceDiri;
+//  std::vector<dirichletBC> faceDiri;
  public:
   elasticitySolver(int tag) : _tag(tag) {}
   void addNodalForces (int iNode, const SVector3 &f)

Solver/quadratureRules.h

+//
+// C++ Interface: quadratureRules
+//
+// Description: 
+//
+//
+// Author:  <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+#ifndef _QUADRATURERULES_H_
+#define _QUADRATURERULES_H_
+
+#include "Gauss.h"
+#include "MElement.h"
+
+class QuadratureBase
+{
+  public : 
+  virtual ~QuadratureBase(){}
+  virtual int getIntPoints(MElement *e,IntPt **GP) =0;
+};
+
+
+class GaussQuadrature : public QuadratureBase
+{
+ public : 
+  enum IntegCases {Other,Val,Grad,ValVal,GradGrad};
+ private :
+  int order;
+  IntegCases info;
+ public : 
+  GaussQuadrature(int order_=0):order(order_),info(Other) {}
+  GaussQuadrature(IntegCases info_):order(0),info(info_) {}
+  virtual ~GaussQuadrature(){}
+  int getIntPoints(MElement *e,IntPt **GP)
+  {
+    int integrationOrder;
+    int npts;
+    int geoorder=e->getPolynomialOrder();    
+    switch(info)
+    {
+    case Other :
+      integrationOrder = order;
+      break;
+    case Val :
+      integrationOrder=geoorder+1;
+      break;
+    case Grad :
+      integrationOrder=geoorder;
+      break;
+    case ValVal :
+      integrationOrder=2*geoorder;
+      break;
+    case GradGrad :
+      integrationOrder=3*(geoorder-1);
+      break;
+    default : integrationOrder=1;
+    }
+    e->getIntegrationPoints(integrationOrder, &npts, GP);
+    return npts;
+  }
+};
+
+
+#endif //_QUADRATURERULES_H_

Solver/solverAlgorithms.h

@@ -15,71 +15,140 @@
 #define _SOLVERALGORITHMS_H_
 
 
-
+#include "dofManager.h"
 #include "terms.h"
+#include "quadratureRules.h"
+#include "MVertex.h"
+#include <iostream>
 
-
-
-template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space1,Iterator itbegin,Iterator itend,Assembler &assembler) // symmetric
+template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler) // symmetric
 {
-    fullMatrix<double> localMatrix;
+  fullMatrix<typename Assembler::dataMat> localMatrix;
   std::vector<Dof> R;
   for (Iterator it = itbegin;it!=itend; ++it)
   {
     MElement *e = *it;
     R.clear();
-      int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ; // todo: to be given from outside
-      int npts=0;
     IntPt *GP;
-      e->getIntegrationPoints(integrationOrder, &npts, &GP);
+    int npts=integrator.getIntPoints(e,&GP);
     term.get(e,npts,GP,localMatrix);
-      space1.getKeys(e,R);
+    space.getKeys(e,R);
     assembler.assemble(R, localMatrix);
   }
 }
 
-template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space1,MElement *e,Assembler &assembler) // symmetric
+template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler) // symmetric
 {
-    fullMatrix<double> localMatrix;
+  fullMatrix<typename Assembler::dataMat> localMatrix;
   std::vector<Dof> R;
-    int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ; // todo: to be given from outside
-    int npts=0;
   IntPt *GP;
-    e->getIntegrationPoints(integrationOrder, &npts, &GP);
+  int npts=integrator.getIntPoints(e,&GP);
   term.get(e,npts,GP,localMatrix);
-    space1.getKeys(e,R);
+  space.getKeys(e,R);
   assembler.assemble(R, localMatrix);
 }
 
-template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space1,Iterator itbegin,Iterator itend,Assembler &assembler)
+template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler)
 {
-    fullVector<double> localVector;
+  fullVector<typename Assembler::dataMat> localVector;
   std::vector<Dof> R;
   for (Iterator it = itbegin;it!=itend; ++it)
   {
     MElement *e = *it;
     R.clear();
-      int integrationOrder = 2* (e->getPolynomialOrder()) ; // todo: to be given from outside
-      int npts=0;
     IntPt *GP;
-      e->getIntegrationPoints(integrationOrder, &npts, &GP);
+    int npts=integrator.getIntPoints(e,&GP);
     term.get(e,npts,GP,localVector);
-      space1.getKeys(e,R);
+    space.getKeys(e,R);
     assembler.assemble(R, localVector);
   }
 }
 
-template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space1,MElement *e,Assembler &assembler)
+template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler)
 {
-    fullVector<double> localVector;
+  fullVector<typename Assembler::dataMat> localVector;
   std::vector<Dof> R;
-    int integrationOrder = 2* (e->getPolynomialOrder()) ; // todo: to be given from outside
-    int npts=0;
   IntPt *GP;
-    e->getIntegrationPoints(integrationOrder, &npts, &GP);
+  int npts=integrator.getIntPoints(e,&GP);
   term.get(e,npts,GP,localVector);
-    space1.getKeys(e,R);
+  space.getKeys(e,R);
   assembler.assemble(R, localVector);
 }
 
+
+template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &dofs,std::vector<typename Assembler::dataVec> &vals)
+{
+  int nbff=dofs.size();
+  for (int i=0;i<nbff;++i)
+  {
+    assembler.fixDof(dofs[i],vals[i]);
+  }
+}
+
+class FilterDof
+{
+ public:
+  virtual bool operator()(Dof key) {return true;}
+};
+
+class FilterDofComponent :public FilterDof
+{
+  int comp;
+ public :
+  FilterDofComponent(int comp_):comp(comp_) {}
+  virtual bool operator()(Dof key)
+  {
+    int type=key.getType();
+    int icomp,iphys;
+    Dof::getTwoIntsFromType(type, iphys, icomp);
+    if (icomp==comp) return true;
+    return false;
+  }
+};
+
+
+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)
+  {
+    MElement *e=*it;
+    std::vector<MVertex*> tabV;
+    int nv=e->getNumVertices();
+    std::vector<Dof> R;
+    space.getKeys(e,R);
+    tabV.reserve(nv);
+    for (int i=0;i<nv;++i) tabV.push_back(e->getVertex(i));
+    
+    for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
+    {
+      Dof key=*itd;
+      if (filter(key))
+      {
+        for (int i=0;i<nv;++i)
+        {
+          if (tabV[i]->getNum()==key.getEntity())
+          {
+            assembler.fixDof(key, fct(tabV[i]->x(),tabV[i]->y(),tabV[i]->z()));
+            break;
+          }
+        }
+      }
+    }
+  }
+}
+
+
+template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
+{
+ for (Iterator it=itbegin;it!=itend;++it)
+  {
+    MElement *e=*it;
+    std::vector<Dof> R;
+    space.getKeys(e,R);
+    int nbdofs=R.size();
+    for (int i=0;i<nbdofs;++i) assembler.numberDof(R[i]);
+  }
+}
+
+
 #endif// _SOLVERALGORITHMS_H_

Solver/terms.h

@@ -19,6 +19,7 @@
 #include <iterator>
 #include "Numeric.h"
 #include "functionSpace.h"
+#include "groupOfElements.h"
 
 
 // evaluation of a field ???
@@ -259,5 +260,32 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
   }
 };
 
+template<class S1> class LoadTermSF : public LinearTerm<S1>
+{
+  simpleFunction<typename S1::ValType> Load;
+ public : 
+  LoadTermSF(S1& space1_,simpleFunction<typename S1::ValType> Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
+  virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
+  {
+    double nbFF=LinearTerm<S1>::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 S1::ValType> Vals;
+      LinearTerm<S1>::space1.f(ele,u, v, w, Vals);
+      for (int j = 0; j < nbFF ; ++j)
+      {
+        m(j)+=dot(Vals[j],Load(u,v,w))*weight*detJ;
+      }
+    }
+  }
+};
+
+
+
 
 #endif// _TERMS_H_