diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9fda9a6f626f0384078b4709770dfe1f5affd9a7..5bebeb67592241af5792190b0c9ee252b4f6bac3 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/Numeric/gmshLinearSystemCSR.h b/Numeric/gmshLinearSystemCSR.h
index 63d536dd0e19e2a1500dedf18c7e47b191aef0ec..735d6139aa71e5d4916b994e9ad6b7e19d31fe44 100644
--- a/Numeric/gmshLinearSystemCSR.h
+++ b/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;
diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h
index 10b1326521d9f455801e379babf7e86143bd10ce..afffe4ce9c85f1ebaa231b9e46d6442ca63fdc6d 100644
--- a/Numeric/gmshLinearSystemGmm.h
+++ b/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);
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index f903ef6583175287b1b2a16882fed8b29df7d8a0..cc7b83f2f4805df8d3a43fd8bd39617d714cb01d 100644
--- a/Solver/CMakeLists.txt
+++ b/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}")
diff --git a/Solver/dofManager.cpp b/Solver/dofManager.cpp
deleted file mode 100644
index 61ee97b1e99354bd340dd269379f601a9c8296e5..0000000000000000000000000000000000000000
--- a/Solver/dofManager.cpp
+++ /dev/null
@@ -1,24 +0,0 @@
-#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;
-}
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 51e6bae57cb7468a570515e27cd5133e7c1f9724..f39bf62ed78cad49eb2d668b42430d40e5e9b777 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -9,74 +9,165 @@
#include <vector>
#include <string>
#include <map>
-
-class gmshLinearSystem;
-
-class Dof{
- private:
- // v(x) = \sum_f \sum_i v_{fi} s^f_i(x)
- long int _entity; // "i": node, edge, group, etc.
- int _type; // "f": basis function type index, etc.
-
- public:
+#include "MVertex.h"
+#include "linearSystem.h"
+
+namespace gsolver {
+
+ class Dof{
+ private:
+ // v(x) = \sum_f \sum_i v_{fi} s^f_i(x)
+ long int _entity; // "i": node, edge, group, etc.
+ int _type; // "f": basis function type index, etc.
+
+ public:
Dof(long int entity, int type) : _entity(entity), _type(type) {}
- inline long int getEntity() const { return _entity; }
- inline int getType() const { return _type; }
-};
-
-template<class dataVec, class dataMat>
-class DofAffineConstraint{
- public:
- std::vector<std::pair<Dof, dataMat> > linear;
- dataVec shift;
-};
-
-// data=float, double, complex<double>, smallVec<double>, smallVec<complex<double> >...
-template <class dataVec, class dataMat>
-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;
-
- // fixations on full blocks, treated by eliminating equations:
- // DofVec = dataVec
- std::map<Dof, dataVec> fixed, initial;
-
- // numbering of unknown dof blocks
- std::map<Dof, long int> unknown;
-
- // 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);
-};
-
-/*
-class termOfFormulation{
- public:
- // hard-coded OR not!
-};
-
-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
-
+ 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:
+ std::vector<std::pair<Dof, dataMat> > linear;
+ dataVec shift;
+ };
+
+ // data=float, double, complex<double>, smallVec<double>, smallVec<complex<double> >...
+ template <class dataVec, class dataMat>
+ 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> constraints;
+
+ // fixations on full blocks, treated by eliminating equations:
+ // DofVec = dataVec
+ std::map<Dof, dataVec> fixed, initial;
+
+ // numbering of unknown dof blocks
+ std::map<Dof, long int> unknown;
+
+ // associatations
+ std::map<Dof, Dof> associatedWith;
+
+ // linearSystems
+ std::map<std::string, linearSystem<dataMat>* > _linearSystems;
+ linearSystem<dataMat>* _current;
+
+ private:
+ public:
+ 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();
+ }
+ };
+}
+#endif
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1545d4cf00d88348b97be30040de3157ce176d99
--- /dev/null
+++ b/Solver/elasticityTerm.cpp
@@ -0,0 +1,104 @@
+// 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;
+ }
+ }
+ }
+}
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6a9fbc3f38d4424713454d1d7dcb9b01010e99e
--- /dev/null
+++ b/Solver/elasticityTerm.h
@@ -0,0 +1,43 @@
+// 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
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..e274702104443a949fe896bcfffb958c62fa7d98
--- /dev/null
+++ b/Solver/linearSystem.h
@@ -0,0 +1,30 @@
+// 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
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5817c1aacffa687112ae16de3e25c3d65a2708cc
--- /dev/null
+++ b/Solver/linearSystemCSR.cpp
@@ -0,0 +1,316 @@
+// 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
+
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
new file mode 100644
index 0000000000000000000000000000000000000000..1dcaae4e714292fe419aa299ef9477afedfe487a
--- /dev/null
+++ b/Solver/linearSystemCSR.h
@@ -0,0 +1,144 @@
+// 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
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
new file mode 100644
index 0000000000000000000000000000000000000000..85ab46429434241497842b5a14839346a4532114
--- /dev/null
+++ b/Solver/linearSystemFull.h
@@ -0,0 +1,85 @@
+// 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
diff --git a/Solver/linearSystemGMM.h b/Solver/linearSystemGMM.h
new file mode 100644
index 0000000000000000000000000000000000000000..822d28fc67079d164089b8a8a9d71a465c818076
--- /dev/null
+++ b/Solver/linearSystemGMM.h
@@ -0,0 +1,122 @@
+// 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
diff --git a/Solver/linearSystemTAUCS.cpp b/Solver/linearSystemTAUCS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4101ca97900c0a8b0561726a51a349259f2d4c2a
--- /dev/null
+++ b/Solver/linearSystemTAUCS.cpp
@@ -0,0 +1,62 @@
+// 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
diff --git a/Solver/linearSystemTAUCS.h b/Solver/linearSystemTAUCS.h
new file mode 100644
index 0000000000000000000000000000000000000000..73a43686bd1e5e382f797a664ce63fced449cce8
--- /dev/null
+++ b/Solver/linearSystemTAUCS.h
@@ -0,0 +1,34 @@
+// 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
diff --git a/Solver/termOfFormulation.h b/Solver/termOfFormulation.h
new file mode 100644
index 0000000000000000000000000000000000000000..7eaa77d9593762cbeff3ee19ebcaa348e852afbe
--- /dev/null
+++ b/Solver/termOfFormulation.h
@@ -0,0 +1,137 @@
+// 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