diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 09ff0d3b1d375e6389b645a3143f1b73a3eef2c6..08e1c2df2b919242706ed734802effa2235689d8 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -199,11 +199,11 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
Msg::Info("Parametrizing Surface %d coordinate %d (ITER %d)", tag(), _isU, ITER);
#ifdef HAVE_GMM
- gmshLinearSystemGmm lsys;
+ gmshLinearSystemGmm<double> lsys;
lsys.setPrec(1.e-10);
- //lsys.setNoisy(2);
+ if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
#else
- gmshLinearSystemFull lsys;
+ gmshLinearSystemFull<double> lsys;
#endif
gmshAssembler myAssembler(&lsys);
gmshGradientBasedDiffusivity diffusivity(coordinates);
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index a26ff2ba29e2c583831e3a019baf3407ea0b76e2..09a4e5557e05e6cd2ea050489653fe73edf16fb6 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -94,7 +94,7 @@ void gmshHighOrderSmoother::smooth(GRegion *gr)
void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
{
- gmshLinearSystemGmm *lsys = new gmshLinearSystemGmm;
+ gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
gmshAssembler myAssembler(lsys);
gmshElasticityTerm El(0,1.0,.333,getTag());
@@ -433,7 +433,7 @@ void localHarmonicMapping(GModel *gm,
// std::vector<SPoint2> &ep4
std::vector<MVertex*> &e) {
- gmshLinearSystemGmm *lsys = new gmshLinearSystemGmm;
+ gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
gmshAssembler myAssembler(lsys);
gmshFunction f(1.0);
gmshLaplaceTerm Laplace (gm,&f,0);
@@ -448,7 +448,7 @@ void localHarmonicMapping(GModel *gm,
Laplace.addToMatrix(myAssembler,t2);
lsys->systemSolve();
- gmshLinearSystemGmm *lsys1 = new gmshLinearSystemGmm;
+ gmshLinearSystemGmm<double> *lsys1 = new gmshLinearSystemGmm<double>;
gmshAssembler myAssembler1(lsys1);
gmshLaplaceTerm Laplace1 (gm,&f,1);
diff --git a/Numeric/gmshAssembler.cpp b/Numeric/gmshAssembler.cpp
index f28f2ec96b450c1a6bbffe4bdf8bd3d3fc7ff305..19160913a596db5e9904c35b7deffe07c681ae78 100644
--- a/Numeric/gmshAssembler.cpp
+++ b/Numeric/gmshAssembler.cpp
@@ -5,7 +5,6 @@
#include "MVertex.h"
#include "gmshAssembler.h"
-#include "gmshLinearSystem.h"
void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
MVertex *vC, int iCompC, int iFieldC,
@@ -15,26 +14,29 @@ void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
lsys->allocate(numbering.size());
}
- std::map<gmshDofKey, int>::iterator itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
+ std::map<gmshDofKey, int>::iterator
+ itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
if (itR != numbering.end()){
- std::map<gmshDofKey, int>::iterator itC = numbering.find(gmshDofKey(vC, iCompC, iFieldC));
+ std::map<gmshDofKey, int>::iterator
+ itC = numbering.find(gmshDofKey(vC, iCompC, iFieldC));
if (itC != numbering.end()){
lsys->addToMatrix(itR->second, itC->second, val);
}
else {
- std::map<gmshDofKey, double>::iterator itF = fixed.find(gmshDofKey(vC, iCompC, iFieldC));
+ std::map<gmshDofKey, double>::iterator
+ itF = fixed.find(gmshDofKey(vC, iCompC, iFieldC));
if (itF != fixed.end()){
lsys->addToRightHandSide(itR->second, -val*itF->second);
}
else{
- std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator itConstrC =
- constraints.find(gmshDofKey(vC, iCompC, iFieldC));
+ std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator
+ itConstrC = constraints.find(gmshDofKey(vC, iCompC, iFieldC));
if (itConstrC != constraints.end()){
for (unsigned int i = 0; i < itConstrC->second.size(); i++){
gmshDofKey &dofKeyConstrC = itConstrC->second[i].first;
double valConstrC = itConstrC->second[i].second;
assemble(vR, iCompR, iFieldR,
- dofKeyConstrC.v,dofKeyConstrC.comp, dofKeyConstrC.field,
+ dofKeyConstrC.v, dofKeyConstrC.comp, dofKeyConstrC.field,
val * valConstrC);
}
}
@@ -42,8 +44,8 @@ void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
}
}
else{
- std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator itConstrR =
- constraints.find(gmshDofKey(vR, iCompR, iFieldR));
+ std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator
+ itConstrR = constraints.find(gmshDofKey(vR, iCompR, iFieldR));
if (itConstrR != constraints.end()){
for (unsigned int i = 0; i < itConstrR->second.size(); i++){
gmshDofKey &dofKeyConstrR = itConstrR->second[i].first;
@@ -60,7 +62,8 @@ void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
double val)
{
if (!lsys->isAllocated())lsys->allocate(numbering.size());
- std::map<gmshDofKey, int>::iterator itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
+ std::map<gmshDofKey, int>::iterator
+ itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
if (itR != numbering.end()){
lsys->addToRightHandSide(itR->second, val);
}
diff --git a/Numeric/gmshAssembler.h b/Numeric/gmshAssembler.h
index 8fac7b745c13a1c3926a7d8b03efd90b3acbb3f1..158e68a7e3dc633ef36cb258d7fb7f4e3fbc4fc7 100644
--- a/Numeric/gmshAssembler.h
+++ b/Numeric/gmshAssembler.h
@@ -8,12 +8,10 @@
#include <map>
#include <vector>
+#include "gmshLinearSystem.h"
class MVertex;
class MElement;
-class gmshLinearSystem;
-
-#include "gmshLinearSystem.h"
struct gmshDofKey{
MVertex *v;
@@ -36,9 +34,9 @@ class gmshAssembler {
std::map<gmshDofKey, int> numbering;
std::map<gmshDofKey, double> fixed;
std::map<gmshDofKey, std::vector<std::pair<gmshDofKey,double> > > constraints;
- gmshLinearSystem *lsys;
+ gmshLinearSystem<double> *lsys;
public:
- gmshAssembler (gmshLinearSystem *l) : lsys(l){}
+ gmshAssembler(gmshLinearSystem<double> *l) : lsys(l) {}
inline void constraintVertex(MVertex*v, int iComp, int iField,
std::vector<MVertex*> &verts,
std::vector<double> &coeffs)
@@ -46,7 +44,7 @@ public:
std::vector<std::pair<gmshDofKey, double> > constraint;
gmshDofKey key(v, iComp, iField);
for (unsigned int i = 0; i < verts.size(); i++){
- gmshDofKey key2 (verts[i],iComp,iField);
+ gmshDofKey key2(verts[i], iComp, iField);
constraint.push_back(std::make_pair(key2, coeffs[i]));
}
constraints[key] = constraint;
@@ -112,8 +110,8 @@ public:
double val);
void assemble(MVertex *vR , int iCompR, int iFieldR,
double val);
- int sizeOfR () const { return numbering.size(); }
- int sizeOfF () const { return fixed.size(); }
+ int sizeOfR() const { return numbering.size(); }
+ int sizeOfF() const { return fixed.size(); }
};
#endif
diff --git a/Numeric/gmshLinearSystem.h b/Numeric/gmshLinearSystem.h
index 4c2eec641c3bfee9fd35bb7589efd4c89e0dfb51..106e4f1f8f5f1d761d0c6755847eada1f0c81f93 100644
--- a/Numeric/gmshLinearSystem.h
+++ b/Numeric/gmshLinearSystem.h
@@ -9,20 +9,20 @@
// A class that encapsulates a linear system solver interface :
// building a sparse matrix, solving a linear system
+template <class scalar>
class gmshLinearSystem {
-public :
+ public :
gmshLinearSystem (){}
- virtual bool isAllocated () const = 0;
- virtual void allocate (int nbRows) = 0;
- virtual ~gmshLinearSystem () {}
- virtual void addToMatrix (int _row, int _col, double val) = 0;
- virtual double getFromMatrix (int _row, int _col) const = 0;
- virtual void addToRightHandSide (int _row, double val) = 0;
- virtual double getFromRightHandSide (int _row) const = 0;
- virtual double getFromSolution (int _row) const = 0;
- virtual void zeroMatrix () = 0;
- virtual void zeroRightHandSide () = 0;
- virtual int systemSolve () = 0;
+ virtual bool isAllocated() const = 0;
+ virtual void allocate(int nbRows) = 0;
+ virtual ~gmshLinearSystem() {}
+ 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/Numeric/gmshLinearSystemFull.h b/Numeric/gmshLinearSystemFull.h
index 054632dcd6d304ac6dcb685f84f5338f9d101691..4ee0ecfa3269b4490966f74afa5fd20cc6d62522 100644
--- a/Numeric/gmshLinearSystemFull.h
+++ b/Numeric/gmshLinearSystemFull.h
@@ -10,61 +10,62 @@
#include "GmshMessage.h"
#include "gmshLinearSystem.h"
-
#include "GmshMatrix.h"
-class gmshLinearSystemFull : public gmshLinearSystem {
- gmshMatrix<double> *_a;
- gmshVector<double> *_b, *_x;
-public :
- gmshLinearSystemFull () : _a(0), _b(0), _x(0){}
- virtual bool isAllocated () const {return _a != 0;}
- virtual void allocate (int _nbRows)
+template <class scalar>
+class gmshLinearSystemFull : public gmshLinearSystem<scalar> {
+ private:
+ gmshMatrix<scalar> *_a;
+ gmshVector<scalar> *_b, *_x;
+ public :
+ gmshLinearSystemFull() : _a(0), _b(0), _x(0){}
+ virtual bool isAllocated() const { return _a != 0; }
+ virtual void allocate(int _nbRows)
{
- if (_a) delete _a;
- if (_x) delete _x;
- if (_b) delete _b;
- _a = new gmshMatrix<double>(_nbRows,_nbRows);
- _b = new gmshVector<double>(_nbRows);
- _x = new gmshVector<double>(_nbRows);
+ if(_a) delete _a;
+ if(_x) delete _x;
+ if(_b) delete _b;
+ _a = new gmshMatrix<scalar>(_nbRows, _nbRows);
+ _b = new gmshVector<scalar>(_nbRows);
+ _x = new gmshVector<scalar>(_nbRows);
}
- virtual ~gmshLinearSystemFull ()
+ virtual ~gmshLinearSystemFull()
{
delete _a;
delete _b;
delete _x;
}
- virtual void addToMatrix (int _row, int _col, double _val)
+ virtual void addToMatrix(int _row, int _col, scalar _val)
{
- if (_val != 0.0) (*_a)(_row, _col) += _val;
+ if(_val != 0.0) (*_a)(_row, _col) += _val;
}
- virtual double getFromMatrix (int _row, int _col) const
+ virtual scalar getFromMatrix(int _row, int _col) const
{
return (*_a)(_row, _col);
}
- virtual void addToRightHandSide (int _row, double _val)
+ virtual void addToRightHandSide(int _row, scalar _val)
{
- if (_val != 0.0) (*_b)(_row)+=_val;
+ if(_val != 0.0) (*_b)(_row) += _val;
}
- virtual double getFromRightHandSide (int _row) const
+ virtual scalar getFromRightHandSide(int _row) const
{
return (*_b)(_row);
}
- virtual double getFromSolution (int _row) const
+ virtual scalar getFromSolution(int _row) const
{
return (*_x)(_row);
}
- virtual void zeroMatrix ()
+ virtual void zeroMatrix()
{
- _a->set_all(0.0);
+ _a->set_all(0.);
}
- virtual void zeroRightHandSide ()
+ virtual void zeroRightHandSide()
{
- for (int i = 0; i < _b->size(); i++) (*_b)(i) = 0;
+ for(int i = 0; i < _b->size(); i++) (*_b)(i) = 0.;
}
- virtual int systemSolve ()
+ virtual int systemSolve()
{
- _a->lu_solve(*_b,*_x);
+ _a->lu_solve(*_b, *_x);
return 1;
}
};
diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h
index 6b5dd077b4f4764a4981abae9ebb23ef4ae80237..461c28d8e218e203cbb00668b1fb9baac791471e 100644
--- a/Numeric/gmshLinearSystemGmm.h
+++ b/Numeric/gmshLinearSystemGmm.h
@@ -16,68 +16,70 @@
#include <gmm.h>
-class gmshLinearSystemGmm : public gmshLinearSystem {
- gmm::row_matrix<gmm::wsvector<double> > *_a;
- std::vector<double> *_b, *_x;
+template <class scalar>
+class gmshLinearSystemGmm : public gmshLinearSystem<scalar> {
+ private:
+ gmm::row_matrix<gmm::wsvector<scalar> > *_a;
+ std::vector<scalar> *_b, *_x;
double _prec;
int _noisy;
int _gmres;
-public :
- gmshLinearSystemGmm () : _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)
+ public:
+ gmshLinearSystemGmm() : _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)
{
- if (_a) delete _a;
- if (_x) delete _x;
- if (_b) delete _b;
- _a = new gmm::row_matrix< gmm::wsvector<double> >(_nbRows,_nbRows);
- _b = new std::vector<double>(_nbRows);
- _x = new std::vector<double>(_nbRows);
+ if(_a) delete _a;
+ if(_x) delete _x;
+ if(_b) delete _b;
+ _a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
+ _b = new std::vector<scalar>(_nbRows);
+ _x = new std::vector<scalar>(_nbRows);
}
- virtual ~gmshLinearSystemGmm ()
+ virtual ~gmshLinearSystemGmm()
{
delete _a;
delete _b;
delete _x;
}
- virtual void addToMatrix (int _row, int _col, double _val)
+ virtual void addToMatrix(int _row, int _col, scalar _val)
{
- if (_val != 0.0) (*_a)(_row, _col) += _val;
+ if(_val != 0.0) (*_a)(_row, _col) += _val;
}
- virtual double getFromMatrix (int _row, int _col) const
+ virtual scalar getFromMatrix (int _row, int _col) const
{
return (*_a)(_row, _col);
}
- virtual void addToRightHandSide (int _row, double _val)
+ virtual void addToRightHandSide(int _row, scalar _val)
{
- if (_val != 0.0) (*_b)[_row]+=_val;
+ if(_val != 0.0) (*_b)[_row] += _val;
}
- virtual double getFromRightHandSide (int _row) const
+ virtual scalar getFromRightHandSide(int _row) const
{
return (*_b)[_row];
}
- virtual double getFromSolution (int _row) const
+ virtual scalar getFromSolution(int _row) const
{
return (*_x)[_row];
}
- virtual void zeroMatrix ()
+ virtual void zeroMatrix()
{
gmm::clear(*_a);
}
- virtual void zeroRightHandSide ()
+ virtual void zeroRightHandSide()
{
- for (unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0;
+ 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 ()
+ 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::rsvector<double> > > P(*_a, 10,0.);
- gmm::ildltt_precond< gmm::row_matrix< gmm::wsvector<double> > > P(*_a, 2, 1.e-10);
+ // gmm::ilutp_precond<gmm::row_matrix<gmm::rsvector<scalar> > > P(*_a, 10,0.);
+ gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 2, 1.e-10);
gmm::iteration iter(_prec); // defines an iteration object, with a max residu of 1E-8
iter.set_noisy(_noisy);
- if (_gmres)gmm::gmres(*_a, *_x, *_b, P, 100, iter); // execute the GMRES algorithm
+ if(_gmres) gmm::gmres(*_a, *_x, *_b, P, 100, iter); // execute the GMRES algorithm
else gmm::cg(*_a, *_x, *_b, P, iter); // execute the CG algorithm
return 1;
}
@@ -85,22 +87,23 @@ public :
#else
+template <class scalar>
class gmshLinearSystemGmm : public gmshLinearSystem {
public :
- gmshLinearSystemGmm ()
+ gmshLinearSystemGmm()
{
- Msg::Error("Gmm++ is not available on this version of gmsh");
+ 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, double val) {}
- virtual double getFromMatrix (int _row, int _col) const { return 0.; }
- virtual void addToRightHandSide (int _row, double val) {}
- virtual double getFromRightHandSide (int _row) const { return 0.; }
- virtual double getFromSolution (int _row) const { return 0.; }
- virtual void zeroMatrix () {}
- virtual void zeroRightHandSide () {}
- virtual int systemSolve () { return 0; }
+ 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; }
};
#endif
diff --git a/Numeric/gmshTermOfFormulation.cpp b/Numeric/gmshTermOfFormulation.cpp
index 554b6b1b48157f68ba51c89e9a5d86aa487434b5..4afb8d88ae8e509fdc7a49d293c234f1d0178092 100644
--- a/Numeric/gmshTermOfFormulation.cpp
+++ b/Numeric/gmshTermOfFormulation.cpp
@@ -9,7 +9,6 @@
#include "GmshMatrix.h"
#include "gmshTermOfFormulation.h"
#include "gmshFunction.h"
-#include "gmshLinearSystem.h"
#include "gmshAssembler.h"
gmshNodalFemTerm::~gmshNodalFemTerm ()
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index ffd3cee46af07b836e0e147d90a9cf3e04f8ef20..966a03896500b916e64103a7fc75d165aac1a8c0 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -6,24 +6,22 @@
#ifndef _GMSH_TERM_OF_FORMULATION_H_
#define _GMSH_TERM_OF_FORMULATION_H_
-class GModel;
-class GEntity;
-class MElement;
-class MVertex;
-class gmshLinearSystem;
-class gmshFunction;
-
#include <math.h>
#include <map>
#include <vector>
#include "GmshMatrix.h"
+class GModel;
+class GEntity;
+class MElement;
+class MVertex;
+class gmshFunction;
class gmshAssembler;
class gmshTermOfFormulation {
-protected:
+ protected:
GModel *_gm;
-public:
+ public:
gmshTermOfFormulation(GModel *gm) : _gm(gm) {}
virtual ~gmshTermOfFormulation(){}
virtual void addToMatrix(gmshAssembler&) const = 0;