*** empty log message ***

CMakeLists.txt

@@ -38,6 +38,7 @@ option(ENABLE_OSMESA "Use OSMesa for offscreen rendering" OFF)
 option(ENABLE_PARSER "Build the GEO file parser" ON)
 option(ENABLE_POST "Build the post-processing module" ON)
 option(ENABLE_QT "Build QT GUI" OFF)
+option(ENABLE_SOLVER "Build Gmsh with the gsolver finite element library" OFF)
 option(ENABLE_TAUCS "Enable Taucs linear algebra solver" ON)
 option(ENABLE_TETGEN "Enable Tetgen mesh generator" ON)
 option(ENABLE_TETGEN_NEW "Enable experimental version of Tetgen" OFF)
@@ -59,7 +60,7 @@ set(GMSH_API
     Numeric/gmshAssembler.h Numeric/gmshTermOfFormulation.h 
     Numeric/gmshLaplace.h Numeric/gmshElasticity.h Numeric/gmshLinearSystem.h
     Numeric/gmshLinearSystemGmm.h Numeric/gmshLinearSystemFull.h 
-    Numeric/gmshFunction.h
+    Numeric/gmshLinearSystemCSR.h Numeric/gmshFunction.h
   Geo/GModel.h Geo/GEntity.h Geo/GPoint.h Geo/GVertex.h Geo/GEdge.h 
     Geo/GFace.h Geo/GRegion.h Geo/GEdgeLoop.h Geo/GEdgeCompound.h 
     Geo/GFaceCompound.h Geo/GRegionCompound.h Geo/MVertex.h Geo/MEdge.h 
@@ -483,6 +484,12 @@ if(ENABLE_TAUCS)
   endif(TAUCS_LIB)
 endif(ENABLE_TAUCS)
 
+if(ENABLE_SOLVER)
+  add_subdirectory(Solver)
+  set(HAVE_GSOLVER TRUE)
+  list(APPEND CONFIG_OPTIONS "gsolver")
+endif(ENABLE_NETGEN)
+
 if(ENABLE_OCC)
   if(WIN32)
     set(OCC_SYS_NAME win32)

Numeric/gmshLinearSystemCSR.h

@@ -45,7 +45,7 @@ class gmshLinearSystemCSR : public gmshLinearSystem<scalar> {
   {
     allocate(0);
   }
-  void addToMatrix ( int il, int ic, double val) 
+  virtual void addToMatrix ( int il, int ic, double val) 
   {
     //    if (sorted)throw;
 

Numeric/gmshLinearSystemGmm.h

@@ -81,8 +81,8 @@ class gmshLinearSystemGmm : public gmshLinearSystem<scalar> {
   void setGmres(int n){ _gmres = n; }
   virtual int systemSolve() 
   {
-    //gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 15, 0.);
-    gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 10, 1.e-10);
+    //gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 25, 0.);
+    gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 30, 1.e-10);
     gmm::iteration iter(_prec);
     iter.set_noisy(_noisy);
     if(_gmres) gmm::gmres(*_a, *_x, *_b, P, 100, iter);

Solver/CMakeLists.txt

@@ -4,7 +4,9 @@
 # bugs and problems to <gmsh@geuz.org>.
 
 set(SRC
-  dofManager.cpp
+  linearSystemCSR.cpp
+  linearSystemTAUCS.cpp
+  elasticityTerm.cpp
 )
 
 append_gmsh_src(Solver "${SRC}")

Solver/dofManager.cpp

-#include "dofManager.h"
-
-template <>
-double DofManager<double, double>::getDofValue(long int entity, int type) const
-{
-  return 1.2345;
-  /*
-     if("find in fixed")
-       return ...
-     else if("find in init")
-       return ...
-     else if("find in constraint")
-       ...
-     else
-       ...
-  */
-}
-
-template <>
-void DofManager<double, double>::createFixedDof(long int entity, int type, 
-                                                const double &value)
-{
-  fixed[Dof(entity, type)] = value;
-}

Solver/dofManager.h

@@ -9,8 +9,10 @@
 #include <vector>
 #include <string>
 #include <map>
+#include "MVertex.h"
+#include "linearSystem.h"
 
-class gmshLinearSystem;
+namespace gsolver {
 
   class Dof{
   private:
@@ -22,8 +24,22 @@ 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 ){
+      return i1 + 10000 * i2;
+    }
+    static void getTwoIntsFromType( int t, int &i1, int &i2){
+      i1 = t % 10000;
+      i2 = t / 10000;
+    }
+    bool operator < ( const Dof & other){
+      if (_entity < other._entity)return true;
+      if (_entity > other._entity)return false;
+      if (_type < other._type)return true;
+      return false;
+    }
   };
 
+
   template<class dataVec, class dataMat>
     class DofAffineConstraint{
   public:
@@ -33,12 +49,12 @@ class DofAffineConstraint{
   
   // data=float, double, complex<double>, smallVec<double>, smallVec<complex<double> >...
   template <class dataVec, class dataMat>
-class DofManager{
+    class dofManager{
   private:
     // general affine constraint on sub-blocks, treated by adding
     // equations:
     //   dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
-  std::map<std::pair<Dof, dataMat>, DofAffineConstraint> constraint;
+    //std::map<std::pair<Dof, dataMat>, DofAffineConstraint> constraints;
     
     // fixations on full blocks, treated by eliminating equations:
     //   DofVec = dataVec    
@@ -50,33 +66,108 @@ class DofManager{
     // associatations
     std::map<Dof, Dof> associatedWith;
     
- private:
-  // linear system(s):
-  std::vector<gmshLinearSystem<dataMat>*> linearSystems;
-
- public:
-  DofManager(/* linAlgOptions &opt */){}
-  dataVec getDofValue(long int entity, int type) const;
-  void createFixedDof(long int entity, int type, const dataVec &value);
-};
+    // linearSystems
+    std::map<std::string, linearSystem<dataMat>* > _linearSystems;
+    linearSystem<dataMat>* _current;  
     
-/*
-class termOfFormulation{
+  private:
   public:
-  // hard-coded OR not!
+    dofManager(linearSystem<dataMat> *l) : _current(l){
+      _linearSystems.push_back(l);
+    }
+    inline void fixDof(long int ent, int type, const dataVec & value)
+    {
+      fixed [Dof(ent, type)] = value;
+    }
+    inline void fixVertex(MVertex*v, int iComp, int iField, const dataVec & value)
+    {
+      fixDof (v->getNum(), Dof::createTypeWithTwoInts(iComp, iField ),value);
+    }
+    inline void numberDof(long int ent, int type){
+      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);
+      if (it == unknown.end()){
+	unsigned int size = unknown.size();
+	unknown[key] = size;
+      }
+    }
+    inline void numberVertex(MVertex*v, int iComp, int iField){
+      numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp,iField ));
+    }
+    inline dataVec getDofValue(int ent, int type) const
+    {
+      Dof key (ent, type);
+      {
+	typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
+	if (it != fixed.end()) return it->second;
+      }
+      {
+	std::map<Dof, int>::const_iterator it = unknown.find(key);
+	if (it != unknown.end())
+	  return _current->getFromSolution(it->second);
+      }
+      return dataVec(0.0);
+    }
+    inline dataVec getVertexValue(MVertex *v, int iComp, int iField) const
+    {
+      getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp,iField ));
+    }
+    
+    inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
+    {
+      if (!_current->isAllocated()) _current->allocate(unknown.size());
+      
+      std::map<Dof, int>::iterator itR = unknown.find(R);
+      if (itR != unknown.end()){
+	std::map<Dof, int>::iterator itC = unknown.find(C);
+	if (itC != unknown.end()){
+	  _current->addToMatrix(itR->second, itC->second, value);
+	}
+	else{
+	  typename std::map<Dof,  dataVec>::iterator itFixed = fixed.find(C);
+	  if (itFixed != fixed.end()){
+	    _current->addToRightHandSide(itR->second, -value*itFixed->second);
+	  }
+	}
+      }
+    }
+
+    inline void assemble(int entR, int typeR, int entC, int typeC, const dataMat &value)
+    {
+      assemble (Dof(entR,typeR),Dof(entC,typeC), value);
+    }
+    
+    inline void assembleVertex(MVertex *vR, int iCompR, int iFieldR,
+			       MVertex *vC, int iCompC, int iFieldC,
+			       const dataMat &value){
+      assemble (vR->getNum(), Dof::createTypeWithTwoInts(iCompR,iFieldR),
+		vC->getNum(), Dof::createTypeWithTwoInts(iCompC,iFieldC),
+		value);
+    } 
+    
+    inline void assemble(const Dof &R, const dataMat &value)
+    {
+      if (!_current->isAllocated())_current->allocate(unknown.size());
+      std::map<Dof, int>::iterator itR = unknown.find(R);
+      if (itR != unknown.end()){
+	_current->addToRightHandSide(itR->second, value);
+      }
+    }
+    inline void assemble(int entR, int typeR, const dataMat &value)
+    {
+      assemble(Dof(entR,typeR));
+    }
+    int sizeOfR() const { return unknown.size(); }
+    int sizeOfF() const { return fixed.size(); }
+    void systemSolve(){
+      _current->systemSolve();
+    }  
+    void systemClear(){
+      _current->zeroMatrix();
+      _current->zeroRightHandSide();
+    }  
   };
-
-class field{
- private:
-  std::string _name;
-
-  // unique mapping between dofs and names:
-  static std::map<int, std::string> _names;
-  static int typeFromName(std::string name);
-  static std::string nameFromType(int type);
-
-  // size of a dof block
-  static int blockSizeFromType(int type);
-
-  // doit pouvoir etre evalu
-
+}
+#endif

Solver/elasticityTerm.cpp

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "elasticityTerm.h"
+#include "Numeric.h"
+
+namespace gsolver {
+  void elasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+  {
+    int nbNodes = e->getNumVertices();
+    int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
+    int npts;
+    IntPt *GP;
+    double jac[3][3];
+    double invjac[3][3];
+    double Grads[256][3], grads[100][3];
+    e->getIntegrationPoints(integrationOrder, &npts, &GP);
+    
+    m.set_all(0.);
+    
+    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;
+    const double C[6][6] =
+      { {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} };
+    
+    gmshMatrix<double> H(6, 6);
+    gmshMatrix<double> B(6, 3 * nbNodes);
+    gmshMatrix<double> BTH(3 * nbNodes, 6);
+    gmshMatrix<double> BT(3 * nbNodes, 6);
+    for (int i = 0; i < 6; i++)
+      for (int j = 0; j < 6; j++)
+	H(i, j) = C[i][j];
+    
+    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 = e->getJacobian(u, v, w, jac);
+      e->getGradShapeFunctions(u, v, w, grads);
+      inv3x3(jac, invjac);
+      B.set_all(0.);
+      BT.set_all(0.);
+      for (int j = 0; j < nbNodes; j++){
+	Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] + 
+	  invjac[0][2] * grads[j][2];
+	Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] + 
+	  invjac[1][2] * grads[j][2];
+	Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] + 
+	  invjac[2][2] * grads[j][2];
+	
+	BT(j, 0) = B(0, j) = Grads[j][0];
+	BT(j, 3) = B(3, j) = Grads[j][1];
+	BT(j, 4) = B(4, j) = Grads[j][2];
+	
+	BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
+	BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
+	BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
+	
+	BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
+	BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
+	BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
+      }
+      BTH.set_all(0.);
+      BTH.gemm(BT, H); 
+      m.gemm(BTH, B, weight * detJ, 1.);
+    } 
+  }
+  
+  void elasticityTerm::elementVector(MElement *e, gmshVector<double> &m) const{
+    int nbNodes = e->getNumVertices();
+    int integrationOrder = 2 * e->getPolynomialOrder();
+    int npts;
+    IntPt *GP;
+    double jac[3][3];
+    double ff[256];
+    e->getIntegrationPoints(integrationOrder, &npts, &GP);
+    
+    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 = e->getJacobian(u, v, w, jac);
+      e->getShapeFunctions(u, v, w, ff);
+      for (int j = 0; j < nbNodes; j++){
+	m(j)           += _volumeForce.x() *ff[j] * weight * detJ *.5;
+	m(j+nbNodes)   += _volumeForce.y() *ff[j] * weight * detJ *.5;
+	m(j+2*nbNodes) += _volumeForce.z() *ff[j] * weight * detJ *.5;
+      }
+    } 
+  }
+}

Solver/elasticityTerm.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _ELASTICITY_H_
+#define _ELASTICITY_H_
+
+#include "termOfFormulation.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "GmshMatrix.h"
+
+namespace gsolver {
+  class elasticityTerm : public femTerm<double,double> {
+  protected:
+    double _E, _nu;
+    int _iField;
+    SVector3 _volumeForce;
+ public:
+    // element matrix size : 3 dofs per vertex
+    virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
+    virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
+    // order dofs in the local matrix :
+    // dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn 
+    Dof getLocalDofR(MElement *e, int iRow) const 
+    {
+      int iCompR = iRow / e->getNumVertices();
+      int ithLocalVertex = iRow % e->getNumVertices();
+      return Dof (e->getVertex(ithLocalVertex)->getNum(),
+		  Dof::createTypeWithTwoInts(iCompR, _iField));
+    }
+  public:
+    elasticityTerm(GModel *gm, double E, double nu, int iField = 1) : 
+      femTerm<double,double>(gm), _E(E), _nu(nu), _iField(iField){}
+    void setVector(const SVector3 &f) {_volumeForce = f;}
+    void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+    void elementVector(MElement *e, gmshVector<double> &m) const;
+  };
+}
+
+#endif

Solver/linearSystem.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _LINEAR_SYSTEM_H_
+#define _LINEAR_SYSTEM_H_
+
+// A class that encapsulates a linear system solver interface :
+// building a sparse matrix, solving a linear system
+namespace gsolver {
+  template <class scalar>
+    class linearSystem {
+  public :
+    linearSystem (){}
+    virtual ~linearSystem (){}
+    virtual bool isAllocated() const = 0;
+    virtual void allocate(int nbRows) = 0;
+    virtual void clear() = 0;
+    virtual void addToMatrix(int _row, int _col, scalar val) = 0;
+    virtual scalar getFromMatrix(int _row, int _col) const = 0;
+    virtual void addToRightHandSide(int _row, scalar val) = 0;
+    virtual scalar getFromRightHandSide(int _row) const = 0;
+    virtual scalar getFromSolution(int _row) const = 0;
+    virtual void zeroMatrix() = 0;
+    virtual void zeroRightHandSide() = 0;
+    virtual int systemSolve() = 0;
+  };
+}
+#endif

Solver/linearSystemCSR.cpp

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "linearSystemCSR.h"
+
+#define SWAP(a,b)  temp=(a);(a)=(b);(b)=temp;
+#define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
+namespace gsolver {
+  void *CSRMalloc(size_t size)
+  {
+    void *ptr;
+    if (!size) return(NULL);
+    ptr = malloc(size);
+    return(ptr);
+  }
+  
+  void *CSRRealloc(void *ptr, size_t size)
+  {
+    if (!size) return(NULL);
+    ptr = realloc(ptr,size);
+    return(ptr);
+  }
+  
+  void CSRList_Realloc(CSRList_T *liste,int n)
+  {
+    char* temp;
+    if (n <= 0) return;
+    if (liste->array == NULL) {
+      liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
+      liste->array = (char *)CSRMalloc(liste->nmax * liste->size);
+    }
+    else {
+      if (n > liste->nmax) {
+	liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
+	temp = (char *)CSRRealloc(liste->array, liste->nmax * liste->size);
+	liste->array = temp;
+      }
+    }
+  }
+  
+  
+  CSRList_T *CSRList_Create(int n, int incr, int size)
+  {
+    CSRList_T *liste;
+    
+    if (n <= 0)  n = 1 ;
+    if (incr <= 0) incr = 1;
+    
+    liste = (CSRList_T *)CSRMalloc(sizeof(CSRList_T));
+    
+    liste->nmax    = 0;
+    liste->incr    = incr;
+    liste->size    = size;
+    liste->n       = 0;
+    liste->isorder = 0;
+    liste->array   = NULL;
+    
+    CSRList_Realloc(liste,n);
+    return(liste);
+  }
+  
+  void CSRList_Delete(CSRList_T *liste)
+  {
+    if (liste != 0) {
+      free(liste->array);
+      free(liste);
+    }
+  }
+  
+  void CSRList_Add(CSRList_T *liste, void *data)
+  {
+    liste->n++;
+    
+    CSRList_Realloc(liste,liste->n);
+    liste->isorder = 0;
+    memcpy(&liste->array[(liste->n - 1) * liste->size],data,liste->size);
+  }
+  
+  int CSRList_Nbr(CSRList_T *liste)
+  {
+    return(liste->n);
+  }
+  
+  template<>
+  void linearSystemCSR<double>::allocate(int _nbRows)
+  {
+    if(a_) {
+      CSRList_Delete(a_);
+      CSRList_Delete(ai_);
+      CSRList_Delete(ptr_);
+      CSRList_Delete(jptr_);
+      delete _x;
+      delete _b;
+      delete something;
+    }
+    
+    if(_nbRows == 0){
+      a_ = 0; 
+      ai_ = 0; 
+      ptr_ = 0; 
+      jptr_ = 0; 
+      _b = 0;
+      _x = 0;
+      something = 0;
+      return;
+    }
+    
+    a_    = CSRList_Create (_nbRows, _nbRows, sizeof(double));
+    ai_   = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
+    ptr_  = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
+    jptr_ = CSRList_Create (_nbRows+1, _nbRows, sizeof(INDEX_TYPE));
+    
+    something = new char[_nbRows];
+    
+    for (int i=0;i<_nbRows;i++)something[i] = 0;
+    
+    _b = new std::vector<double>(_nbRows);
+    _x = new std::vector<double>(_nbRows);
+  }
+  
+  const int NSTACK   = 50;
+  const unsigned int M_sort2  = 7;
+  
+  static void free_ivector(int *v, long nl, long nh)
+  {
+    // free an int vector allocated with ivector() 
+    free((char*)(v+nl-1));
+  }
+  
+  static int *ivector(long nl, long nh)
+  {
+    // allocate an int vector with subscript range v[nl..nh] 
+    int *v;  
+    v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
+    if (!v) fprintf(stderr, "allocation failure in ivector()\n");
+    return v-nl+1;
+  }
+  
+  static int cmpij(INDEX_TYPE ai,INDEX_TYPE aj,INDEX_TYPE bi,INDEX_TYPE bj)
+  {
+    if(ai<bi)return -1;
+    if(ai>bi)return 1;
+    if(aj<bj)return -1;
+    if(aj>bj)return 1;
+    return 0;
+  }
+  
+  template <class scalar>
+  void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
+  {
+    unsigned long i,ir=n,j,k,l=1;
+    int *istack,jstack=0;
+    INDEX_TYPE tempi;
+    scalar a,temp;
+    int    b,c;
+    
+    istack=ivector(1,NSTACK);
+    for (;;) {
+      if (ir-l < M_sort2) {
+	for (j=l+1;j<=ir;j++) {
+	  a=arr[j -1];
+	  b=ai[j -1];
+	  c=aj[j -1];
+	  for (i=j-1;i>=1;i--) {
+	    if (cmpij(ai[i -1],aj[i -1],b,c) <= 0) break;
+	    arr[i+1 -1]=arr[i -1];
+	    ai[i+1 -1]=ai[i -1];
+	    aj[i+1 -1]=aj[i -1];
+	  }
+	  arr[i+1 -1]=a;
+	  ai[i+1 -1]=b;
+	  aj[i+1 -1]=c;
+	}
+	if (!jstack) {
+	  free_ivector(istack,1,NSTACK);
+	  return;
+	}
+	ir=istack[jstack];
+	l=istack[jstack-1];
+	jstack -= 2;
+      } 
+      else {
+	k=(l+ir) >> 1;
+	SWAP(arr[k -1],arr[l+1 -1])
+	  SWAPI(ai[k -1],ai[l+1 -1])
+	  SWAPI(aj[k -1],aj[l+1 -1])
+	  if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[ir -1],aj[ir -1])>0){
+	    SWAP(arr[l+1 -1],arr[ir -1])
+	      SWAPI(ai[l+1 -1],ai[ir -1])
+	      SWAPI(aj[l+1 -1],aj[ir -1])
+	      }
+	if (cmpij(ai[l -1],aj[l -1],ai[ir -1],aj[ir -1])>0){
+	  SWAP(arr[l -1],arr[ir -1])
+	    SWAPI(ai[l -1],ai[ir -1])
+	    SWAPI(aj[l -1],aj[ir -1])
+	    }
+	if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[l -1],aj[l -1])>0){
+	  SWAP(arr[l+1 -1],arr[l -1])
+	    SWAPI(ai[l+1 -1],ai[l -1])
+	    SWAPI(aj[l+1 -1],aj[l -1])
+	    }
+	i=l+1;
+	j=ir;
+	a=arr[l -1];
+	b=ai[l -1];
+	c=aj[l -1];
+	for (;;) {
+	  do i++; while (cmpij(ai[i -1],aj[i -1],b,c) < 0);
+	  do j--; while (cmpij(ai[j -1],aj[j -1],b,c) > 0);
+	  if (j < i) break;
+	  SWAP(arr[i -1],arr[j -1])
+	    SWAPI(ai[i -1],ai[j -1])
+	    SWAPI(aj[i -1],aj[j -1])
+	    }
+	arr[l -1]=arr[j -1];
+	arr[j -1]=a;
+	ai[l -1]=ai[j -1];
+	ai[j -1]=b;
+	aj[l -1]=aj[j -1];
+	aj[j -1]=c;
+	jstack += 2;
+	if (jstack > NSTACK) {
+	  Msg::Fatal("NSTACK too small while sorting the columns of the matrix");
+	  throw;
+	}
+	if (ir-i+1 >= j-l) {
+	  istack[jstack]=ir;
+	  istack[jstack-1]=i;
+	  ir=j-1;
+	} 
+	else {
+	  istack[jstack]=j-1;
+	  istack[jstack-1]=l;
+	  l=i;
+	}
+      }
+    }
+  }
+  
+  template <class scalar>
+  void sortColumns(int NbLines, 
+		   int nnz, 
+		   INDEX_TYPE *ptr, 
+		   INDEX_TYPE *jptr, 
+		   INDEX_TYPE *ai, 
+		   scalar *a)
+  {
+    // replace pointers by lines
+    int *count = new int [NbLines];
+    
+    for(int i=0;i<NbLines;i++){
+      count[i] = 0;
+      INDEX_TYPE _position =  jptr[i];
+      while(1){
+	count[i]++;
+	INDEX_TYPE _position_temp = _position;
+	_position = ptr[_position];
+	ptr[_position_temp] = i;
+	if (_position == 0) break;
+      }
+    }   
+    _sort2_xkws<double>(nnz,a,ptr,ai);   
+    jptr[0] = 0;
+    for(int i=1;i<=NbLines;i++){
+      jptr[i] = jptr[i-1] + count[i-1];
+    }
+    
+    for(int i=0;i<NbLines;i++){
+      for (int j= jptr[i] ; j<jptr[i+1]-1 ; j++){
+	ptr[j] = j+1;
+      }
+      ptr[jptr[i+1]] = 0;
+    }
+    
+    delete[] count;
+  }
+}
+  
+#if defined(HAVE_GMM)
+#include "gmm.h"
+namespace gsolver {
+  template<>
+  int linearSystemCSRGmm<double>::systemSolve()
+  {
+    if (!sorted)
+      sortColumns(_b->size(),
+		  CSRList_Nbr(a_),
+		  (INDEX_TYPE *) ptr_->array,
+		  (INDEX_TYPE *) jptr_->array, 
+		  (INDEX_TYPE *) ai_->array, 
+		  (double*) a_->array);
+    sorted = true;
+    
+    gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0>  
+      ref((double*)a_->array, (INDEX_TYPE *) ai_->array,
+	  (INDEX_TYPE *)jptr_->array, _b->size(), _b->size());
+    gmm::csr_matrix<double,0> M;
+    M.init_with(ref);
+    
+    gmm::ildltt_precond<gmm::csr_matrix<double, 0> > P(M, 10, 1.e-10);
+    gmm::iteration iter(_prec);
+    iter.set_noisy(_noisy);
+    if(_gmres) gmm::gmres(M, *_x, *_b, P, 100, iter);
+    else gmm::cg(M, *_x, *_b, P, iter);
+    return 1;
+  }
+}
+#endif
+

Solver/linearSystemCSR.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _LINEAR_SYSTEM_CSR_H_
+#define _LINEAR_SYSTEM_CSR_H_
+
+#include <vector>
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "linearSystem.h"
+
+namespace gsolver {  
+  typedef int INDEX_TYPE ;  
+  typedef struct {
+    int nmax;
+    int size;
+    int incr;
+    int n;
+    int isorder;
+    char *array;
+  } CSRList_T;
+  
+  void CSRList_Add(CSRList_T *liste, void *data);
+  int  CSRList_Nbr(CSRList_T *liste);
+  
+  template <class scalar>
+  class linearSystemCSR : public linearSystem<scalar> {
+  protected:
+    bool sorted;
+    char *something;
+    CSRList_T *a_,*ai_,*ptr_,*jptr_; 
+    std::vector<scalar> *_b, *_x;
+  public:
+    linearSystemCSR()
+      : sorted(false), a_(0) {}
+    virtual bool isAllocated() const { return a_ != 0; }
+    virtual void allocate(int) ;
+    virtual void clear()
+    {
+      allocate(0);
+    }
+    virtual ~linearSystemCSR()
+    {
+      allocate(0);
+    }
+    virtual void addToMatrix ( int il, int ic, double val) 
+    {
+      //    if (sorted)throw;
+      
+      INDEX_TYPE  *jptr  = (INDEX_TYPE*) jptr_->array;
+      INDEX_TYPE  *ptr   = (INDEX_TYPE*) ptr_->array;
+      INDEX_TYPE  *ai    = (INDEX_TYPE*) ai_->array;
+      scalar      *a     = ( scalar * ) a_->array;
+      
+      INDEX_TYPE  position_ = jptr[il];
+      
+      if(something[il]) {
+	while(1){
+	  if(ai[position_] == ic){
+	    a[position_] += val;
+	    //      if (il == 0)    printf("FOUND %d %d %d\n",il,ic,position_);
+	    return;
+	  }
+	  if (ptr[position_] == 0)break;
+	  position_ = ptr[position_];
+	}
+      }  
+      
+      INDEX_TYPE zero = 0;
+      CSRList_Add (a_, &val);
+      CSRList_Add (ai_, &ic);
+      CSRList_Add (ptr_, &zero);
+      // The pointers may have been modified
+      // if there has been a reallocation in CSRList_Add  
+      
+      ptr = (INDEX_TYPE*) ptr_->array;
+      ai  = (INDEX_TYPE*) ai_->array;
+      a   = (scalar*) a_->array;
+      
+      INDEX_TYPE n = CSRList_Nbr(a_) - 1;
+      
+      if(!something[il]) {
+	jptr[il] = n;
+	something[il] = 1;      
+      }
+      else ptr[position_] = n;
+    }
+    
+    virtual scalar getFromMatrix (int _row, int _col) const
+    {
+      throw;
+    }
+    virtual void addToRightHandSide(int _row, scalar _val) 
+    {
+      if(_val != 0.0) (*_b)[_row] += _val;
+    }
+    virtual scalar getFromRightHandSide(int _row) const 
+    {
+      return (*_b)[_row];
+    }
+    virtual scalar getFromSolution(int _row) const
+    {
+      return (*_x)[_row];
+    }
+    virtual void zeroMatrix()
+    {
+      int N=CSRList_Nbr(a_);
+      scalar *a = (scalar*) a_->array;
+      for (int i=0;i<N;i++)a[i]=0;
+    }
+    virtual void zeroRightHandSide() 
+    {
+      for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
+    }
+  };
+
+  template <class scalar>
+  class linearSystemCSRGmm : public linearSystemCSR<scalar> {
+  private:
+    double _prec;
+    int _noisy, _gmres;
+  public:
+    linearSystemCSRGmm()
+      : _prec(1.e-8), _noisy(0), _gmres(0) {}
+    virtual ~linearSystemCSRGmm()
+    {}
+    void setPrec(double p){ _prec = p; }
+    void setNoisy(int n){ _noisy = n; }
+    void setGmres(int n){ _gmres = n; }
+    virtual int systemSolve() 
+#if defined(HAVE_GMM)
+      ;
+#else
+    {
+      Msg::Error("Gmm++ is not available in this version of Gmsh");
+      return 0;
+    }
+#endif
+  };
+}
+
+#endif

Solver/linearSystemFull.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _LINEAR_SYSTEM_FULL_H_
+#define _LINEAR_SYSTEM_FULL_H_
+
+// Interface to a linear system with a FULL matrix
+
+#include "GmshMessage.h"
+#include "linearSystem.h"
+#include "GmshMatrix.h"
+#include <stdlib.h>
+#include <set>
+
+namespace gsolver {
+  template <class scalar>
+  class linearSystemFull : public linearSystem<scalar> {
+  private:
+    gmshMatrix<scalar> *_a;
+    gmshVector<scalar> *_b, *_x;
+  public :
+    linearSystemFull() : _a(0), _b(0), _x(0){}
+    virtual bool isAllocated() const { return _a != 0; }
+    virtual void allocate(int _nbRows)
+    {
+      clear();
+      _a = new gmshMatrix<scalar>(_nbRows, _nbRows);
+      _b = new gmshVector<scalar>(_nbRows);
+      _x = new gmshVector<scalar>(_nbRows);
+    }
+    
+    virtual ~linearSystemFull()
+    {
+      clear();
+    }
+    
+    virtual void clear()
+    {
+      if(_a){
+	delete _a;
+	delete _b;
+	delete _x;
+      }
+      _a = 0;
+    }
+    virtual void addToMatrix(int _row, int _col, scalar _val)
+    {
+      if(_val != 0.0) (*_a)(_row, _col) += _val;
+    }
+    virtual scalar getFromMatrix(int _row, int _col) const
+    {
+      return (*_a)(_row, _col);
+    }
+    virtual void addToRightHandSide(int _row, scalar _val)
+    {
+      if(_val != 0.0) (*_b)(_row) += _val;
+    }
+    virtual scalar getFromRightHandSide(int _row) const 
+    {
+      return (*_b)(_row);
+    }
+    virtual scalar getFromSolution(int _row) const 
+    {
+      return (*_x)(_row);
+    }
+    virtual void zeroMatrix() 
+    {
+      _a->set_all(0.);
+    }
+    virtual void zeroRightHandSide()
+    {
+      for(int i = 0; i < _b->size(); i++) (*_b)(i) = 0.;
+    }
+    virtual int systemSolve() 
+    {
+      if (_b->size())
+	_a->lu_solve(*_b, *_x);
+      // _x->print("********* mySol");
+      return 1;
+    }
+  };
+}
+#endif

Solver/linearSystemGMM.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _LINEAR_SYSTEM_GMM_H_
+#define _LINEAR_SYSTEM_GMM_H_
+
+// Interface to GMM++
+
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "linearSystem.h"
+
+#if defined(HAVE_GMM)
+#include <gmm.h>
+namespace gsolver {
+  template <class scalar>
+  class linearSystemGmm : public linearSystem<scalar> {
+  private:
+    gmm::row_matrix<gmm::wsvector<scalar> > *_a;
+    std::vector<scalar> *_b, *_x;
+    double _prec;
+    int _noisy, _gmres;
+  public:
+    linearSystemGmm()
+      : _a(0), _b(0), _x(0), _prec(1.e-8), _noisy(0), _gmres(0) {}
+    virtual bool isAllocated() const { return _a != 0; }
+    virtual void allocate(int _nbRows)
+    {
+      clear();
+      _a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
+      _b = new std::vector<scalar>(_nbRows);
+      _x = new std::vector<scalar>(_nbRows);
+    }
+    virtual ~linearSystemGmm()
+    {
+      clear();
+    }
+    
+    virtual void clear()
+    {
+      if (_a){
+	delete _a;
+	delete _b;
+	delete _x;
+      }
+      _a = 0;
+    }
+    virtual void  addToMatrix(int _row, int _col, scalar _val) 
+    {
+      if(_val != 0.0) (*_a)(_row, _col) += _val;
+    }
+    virtual scalar getFromMatrix (int _row, int _col) const
+    {
+      return (*_a)(_row, _col);
+    }
+    virtual void addToRightHandSide(int _row, scalar _val) 
+    {
+      if(_val != 0.0) (*_b)[_row] += _val;
+    }
+    virtual scalar getFromRightHandSide(int _row) const 
+    {
+      return (*_b)[_row];
+    }
+    virtual scalar getFromSolution(int _row) const
+    {
+      return (*_x)[_row];
+    }
+    virtual void zeroMatrix()
+    {
+      gmm::clear(*_a);
+    }
+    virtual void zeroRightHandSide() 
+    {
+      for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
+    }
+    void setPrec(double p){ _prec = p; }
+    void setNoisy(int n){ _noisy = n; }
+    void setGmres(int n){ _gmres = n; }
+    virtual int systemSolve() 
+    {
+      //gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 25, 0.);
+      gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 30, 1.e-10);
+      gmm::iteration iter(_prec);
+      iter.set_noisy(_noisy);
+      if(_gmres) gmm::gmres(*_a, *_x, *_b, P, 100, iter);
+      else gmm::cg(*_a, *_x, *_b, P, iter);
+      return 1;
+    }
+  };
+}
+
+#else
+
+namespace gsolver {
+  template <class scalar>
+  class linearSystemGmm : public linearSystem<scalar> {
+  public :
+    linearSystemGmm()
+    {
+      Msg::Error("Gmm++ is not available in this version of Gmsh");
+    }
+    virtual bool isAllocated() const { return false; }
+    virtual void allocate(int nbRows) {}
+    virtual void addToMatrix(int _row, int _col, scalar val) {}
+    virtual scalar getFromMatrix(int _row, int _col) const { return 0.; }
+    virtual void addToRightHandSide(int _row, scalar val) {}
+    virtual scalar getFromRightHandSide(int _row) const { return 0.; }
+    virtual scalar getFromSolution(int _row) const { return 0.; }
+    virtual void zeroMatrix() {}
+    virtual void zeroRightHandSide() {}
+    virtual int systemSolve() { return 0; }
+    void setPrec(double p){}
+    virtual void clear(){}
+    void setNoisy(int n){}
+    void setGmres(int n){}
+  }; 
+}
+#endif
+
+#endif

Solver/linearSystemTAUCS.cpp

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "linearSystemTAUCS.h"
+
+#if defined(HAVE_TAUCS)
+extern "C" {
+#include "taucs.h"
+}
+
+namespace gsolver {
+  template <class scalar>
+  void sortColumns(int NbLines, 
+		   int nnz, 
+		   INDEX_TYPE *ptr, 
+		   INDEX_TYPE *jptr, 
+		   INDEX_TYPE *ai, 
+		   scalar *a);
+  
+  template<>
+  int linearSystemCSRTaucs_<double>::systemSolve()
+  {
+    if(!sorted){
+      sortColumns(_b->size(),
+		  CSRList_Nbr(a_),
+		  (INDEX_TYPE *) ptr_->array,
+		  (INDEX_TYPE *) jptr_->array, 
+		  (INDEX_TYPE *) ai_->array, 
+		  (double*) a_->array);
+    }
+    sorted = true;
+    
+    taucs_ccs_matrix myVeryCuteTaucsMatrix;
+    myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m =  _b->size();
+    //myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ptr_->array;
+    //myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)ai_->array;
+    myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ai_->array;
+    myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)jptr_->array;
+    myVeryCuteTaucsMatrix.values.d = (double*) a_->array;
+    myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE; 
+    
+    char* options[] = { "taucs.factor.LLT=true", NULL };  
+    int result = taucs_linsolve(&myVeryCuteTaucsMatrix, 
+				NULL, 
+				1, 
+				&(*_x)[0],
+				&(*_b)[0],
+				options,
+				NULL);                         
+    if (result != TAUCS_SUCCESS){
+      Msg::Error("Taucs Was Not Successfull %d",result);
+    }  
+    return 1;
+  }
+}
+#endif

Solver/linearSystemTAUCS.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _LINEAR_SYSTEM_TAUCS_H_
+#define _LINEAR_SYSTEM_TAUCS_H_
+
+#include "linearSystemCSR.h"
+
+namespace gsolver {
+  template <class scalar>
+    class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
+  private:
+  public:
+    linearSystemCSRTaucs()
+      {}
+    virtual ~linearSystemCSRTaucs()
+      {}
+    virtual void addToMatrix ( int il, int ic, double val) {
+      if (il <= ic)
+	linearSystemCSR<scalar>::addToMatrix(il,ic,val);
+    }
+
+    virtual int systemSolve() 
+#if defined(HAVE_TAUCS)
+;
+#else
+    {throw;}
+#endif
+  };
+}
+
+#endif

Solver/termOfFormulation.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _TERM_OF_FORMULATION_H_
+#define _TERM_OF_FORMULATION_H_
+
+#include <math.h>
+#include <map>
+#include <vector>
+#include "GmshMatrix.h"
+#include "gmshFunction.h"
+#include "dofManager.h"
+#include "GModel.h"
+#include "MElement.h"
+
+namespace gsolver {
+  // a nodal finite element term : variables are always defined at nodes
+  // of the mesh
+  template<class dataVec, class dataMat>
+  class femTerm {
+  protected:
+    GModel *_gm;
+  public:
+    // return the number of columns of the element matrix
+    virtual int sizeOfC(MElement*) const = 0;
+    // return the number of rows of the element matrix
+    virtual int sizeOfR(MElement*) const = 0;
+    // in a given element, return the dof associated to a given row (column)
+    // of the local element matrix
+    virtual Dof getLocalDofR(MElement *e, int iRow) const = 0;
+    // default behavior : symmetric
+    virtual Dof getLocalDofC(MElement *e, int iCol) const
+    {getLocalDofR(e, iCol);}
+  public:
+    femTerm(GModel *gm) : _gm(gm) {}
+    virtual ~femTerm (){}
+    virtual void elementMatrix(MElement *e, gmshMatrix<dataMat> &m) const = 0;
+    virtual void elementVector(MElement *e, gmshVector<dataVec> &m) const {
+      m.scale(0.0);
+    }
+    void addToMatrix(dofManager<dataVec,dataMat> &dm, GEntity *ge) const
+    {
+      for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+	MElement *e = ge->getMeshElement(i);
+	addToMatrix(dm, e);
+      }
+    }
+    void addToMatrix(dofManager<dataVec,dataMat> &dm, MElement *e) const
+    {
+      const int nbR = sizeOfR(e);
+      const int nbC = sizeOfC(e);
+      gmshMatrix<dataMat> localMatrix(nbR, nbC);
+      elementMatrix(e, localMatrix);
+      addToMatrix(dm, localMatrix, e);
+    }
+    void addToMatrix(dofManager<dataVec,dataMat> &dm, 
+		     gmshMatrix<dataMat> &localMatrix, 
+		     MElement *e) const
+    {
+      const int nbR = localMatrix.size1();
+      const int nbC = localMatrix.size2();
+      for (int j = 0; j < nbR; j++){
+	Dof R = getLocalDofR(e, j);
+	for (int k = 0; k < nbC; k++){
+	  Dof C = getLocalDofC(e, k);
+	  dm.assemble(R,C, localMatrix(j, k));
+	}
+      }
+    }
+    void dirichletNodalBC(int physical,
+			  int dim,
+			  int comp,
+			  int field,
+			  const gmshFunction<dataVec> &e,
+			  dofManager<dataVec,dataMat> &dm)
+    {
+      std::vector<MVertex *> v;
+      GModel *m = _gm;
+      m->getMeshVerticesForPhysicalGroup(dim, physical, v);
+      for (unsigned int i = 0; i < v.size(); i++)
+	dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
+    }
+    void neumannNodalBC(int physical,
+			int dim,
+			int comp,
+			int field,
+			const gmshFunction<dataVec> &e,
+			dofManager<dataVec,dataMat> &dm)
+    {
+      std::map<int, std::vector<GEntity*> > groups[4];
+      GModel *m = _gm;
+      m->getPhysicalGroups(groups);
+      std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);  
+      if (it == groups[dim].end()) return;
+      double jac[3][3];
+      double sf[256];
+      for (unsigned int i = 0; i < it->second.size(); ++i){
+	GEntity *ge = it->second[i];
+	for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
+	  MElement *e = ge->getMeshElement(j);
+	  int integrationOrder = 2 * e->getPolynomialOrder();
+	  int nbNodes = e->getNumVertices();
+	  int npts;
+	  IntPt *GP;
+	  e->getIntegrationPoints(integrationOrder, &npts, &GP);  
+	  for (int ip = 0; ip < npts; ip++){
+	    const double u = GP[ip].pt[0];
+	    const double v = GP[ip].pt[1];
+	    const double w = GP[ip].pt[2];
+	    const double weight = GP[ip].weight;
+	    const double detJ = e->getJacobian(u, v, w, jac);
+	    SPoint3 p; e->pnt(u, v, w, p);
+	    e->getShapeFunctions(u, v, w, sf);
+	    const dataVec FCT = fct(p.x(), p.y(), p.z());
+	    for (int k = 0; k < nbNodes; k++){
+	     dm.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
+	    }
+	  }
+	}
+      }
+    }
+    void addToRightHandSide (dofManager<dataVec,dataMat> &dm, GEntity *ge) const {
+      for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+	MElement *e = ge->getMeshElement (i);
+	int nbR = sizeOfR(e);
+	gmshVector<dataVec> V (nbR);
+	elementVector (e, V);
+	// assembly
+	for (int j=0;j<nbR;j++)dm.assemble(getLocalDofR(e,j),V[j]);
+      }
+    }
+  }; // end of class definition
+}// end of namespace
+
+#endif