Skip to content
Snippets Groups Projects
Commit e72c3d5d authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

templatize gmshLinearSystem

parent 5364498f
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 120 additions and 118 deletions
Geo/GFaceCompound.cpp
+ 3
− 3
View file @ e72c3d5d
......@@ -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);
......
Mesh/gmshSmoothHighOrder.cpp
+ 3
− 3
View file @ e72c3d5d
......@@ -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);
......
Numeric/gmshAssembler.cpp
+ 13
− 10
View file @ e72c3d5d
......@@ -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,20 +14,23 @@ 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;
......@@ -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);
}
......
Numeric/gmshAssembler.h
+ 6
− 8
View file @ e72c3d5d
......@@ -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)
......
Numeric/gmshLinearSystem.h
+ 13
− 13
View file @ e72c3d5d
......@@ -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 :
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 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
Numeric/gmshLinearSystemFull.h
+ 29
− 28
View file @ e72c3d5d
......@@ -10,12 +10,13 @@
#include "GmshMessage.h"
#include "gmshLinearSystem.h"
#include "GmshMatrix.h"
class gmshLinearSystemFull : public gmshLinearSystem {
gmshMatrix<double> *_a;
gmshVector<double> *_b, *_x;
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; }
......@@ -24,9 +25,9 @@ public :
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);
_a = new gmshMatrix<scalar>(_nbRows, _nbRows);
_b = new gmshVector<scalar>(_nbRows);
_x = new gmshVector<scalar>(_nbRows);
}
virtual ~gmshLinearSystemFull()
{
......@@ -34,33 +35,33 @@ public :
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;
}
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;
}
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()
{
_a->set_all(0.0);
_a->set_all(0.);
}
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()
{
......
Numeric/gmshLinearSystemGmm.h
+ 46
− 43
View file @ e72c3d5d
......@@ -16,9 +16,11 @@
#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;
......@@ -30,9 +32,9 @@ public :
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);
_a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
_b = new std::vector<scalar>(_nbRows);
_x = new std::vector<scalar>(_nbRows);
}
virtual ~gmshLinearSystemGmm()
{
......@@ -40,23 +42,23 @@ public :
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;
}
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;
}
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];
}
......@@ -66,15 +68,15 @@ public :
}
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()
{
// 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
......@@ -85,19 +87,20 @@ public :
#else
template <class scalar>
class gmshLinearSystemGmm : public gmshLinearSystem {
public :
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 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; }
......
Numeric/gmshTermOfFormulation.cpp
+ 0
− 1
View file @ e72c3d5d
......@@ -9,7 +9,6 @@
#include "GmshMatrix.h"
#include "gmshTermOfFormulation.h"
#include "gmshFunction.h"
#include "gmshLinearSystem.h"
#include "gmshAssembler.h"
gmshNodalFemTerm::~gmshNodalFemTerm ()
......
Numeric/gmshTermOfFormulation.h
+ 7
− 9
View file @ e72c3d5d
......@@ -6,18 +6,16 @@
#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 {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment