diff --git a/Solver/function.cpp b/Solver/function.cpp
index 1f175c6c4e5f9a388d6f640c6739dd8c0d9d7f32..51bd67e50353c795ec097443c1c64c98e8d24484 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -631,3 +631,6 @@ functionConstant *function::getDT() {
_dtFunction = functionConstantNew(0.);
return _dtFunction;
}
+function *getFunctionCoordinates(){
+ return functionCoordinates::get();
+}
diff --git a/Solver/function.h b/Solver/function.h
index b341fcce31fe3408860af32dbf8a845d9bc13417..945b4b11fb3eb3f6f3d4fac862a560012410b28a 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -242,4 +242,5 @@ class functionConstant : public function {
}
void set(double val);
};
+function *getFunctionCoordinates();
#endif
diff --git a/Solver/linearSystem.cpp b/Solver/linearSystem.cpp
index f823902db4734496fa11961746cd017994f04df7..edc45772d0b021b02478299de9208a502e84d804 100644
--- a/Solver/linearSystem.cpp
+++ b/Solver/linearSystem.cpp
@@ -43,3 +43,9 @@ void linearSystem<double>::registerBindings(binding *b){
linearSystemPETScRegisterBindings (b);
#endif
}
+
+void linearSystemBase::setParameter (std::string key, std::string value) {
+ if (isAllocated())
+ Msg::Error("this system is already allocated, parameters cannot be set");
+ _parameters[key] = value;
+}
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index bec6ac2bd4c181f1bfbae3f47abaf1e9c8f94545..687aa14d5ca1d7e37a23c68dcbe5a5b3ea9b3759 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -5,27 +5,35 @@
#ifndef _LINEAR_SYSTEM_H_
#define _LINEAR_SYSTEM_H_
+#include <map>
+#include <string>
// A class that encapsulates a linear system solver interface :
// building a sparse matrix, solving a linear system
class binding;
+class linearSystemBase {
+ protected:
+ std::map<std::string, std::string> _parameters;
+ public :
+ virtual bool isAllocated() const = 0;
+ virtual void allocate(int nbRows) = 0;
+ virtual void clear() = 0;
+ virtual void zeroMatrix() = 0;
+ virtual void zeroRightHandSide() = 0;
+ virtual int systemSolve() = 0;
+ void setParameter (std::string key, std::string value);
+};
template <class scalar>
-class linearSystem {
+class linearSystem : public linearSystemBase {
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, const scalar &val) = 0;
virtual void getFromMatrix(int _row, int _col, scalar &val) const = 0;
virtual void addToRightHandSide(int _row, const scalar &val) = 0;
virtual void getFromRightHandSide(int _row, scalar &val) const = 0;
virtual void getFromSolution(int _row, scalar &val) const = 0;
- virtual void zeroMatrix() = 0;
- virtual void zeroRightHandSide() = 0;
- virtual int systemSolve() = 0;
static void registerBindings (binding*);
virtual double normInfRightHandSide() const = 0;
};
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index cad9f40901fcee78b94c88e798003a9c19038f71..876f9aa74bc115cf7e25a688ffedfd2a1248ebc6 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -3,6 +3,8 @@
#if defined(HAVE_PETSC)
#include "linearSystemPETSc.h"
#include "fullMatrix.h"
+#include <stdlib.h>
+#include "GmshMessage.h"
template <>
void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
@@ -56,6 +58,7 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromRightHandSide(int row,
}
_try(VecRestoreArray(_b, &tmp));
}
+
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullMatrix<PetscScalar> &val) const
{
@@ -75,6 +78,9 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullM
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
{
+ _blockSize = strtol (_parameters["blockSize"].c_str(), NULL, 10);
+ if (_blockSize == 0)
+ Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
clear();
_try(MatCreate(PETSC_COMM_WORLD, &_a));
_try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows * _blockSize, nbRows * _blockSize));
@@ -108,10 +114,9 @@ void linearSystemPETScRegisterBindings(binding *b)
cm->setArgNames(NULL);
cb = b->addClass<linearSystemPETSc<fullMatrix<PetscScalar> > >("linearSystemPETScBlock");
cb->setDescription("A linear system solver, based on PETSc");
- cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> >, int>();
- cm->setDescription ("A new PETScBlock<PetscScalar> solver (we probably should get rid of the blockSize argument)");
+ cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> > >();
+ cm->setDescription ("A new PETScBlock<PetscScalar> solver");
cb->setParentClass<linearSystem<fullMatrix<PetscScalar> > >();
- cm->setArgNames("blockSize", NULL);
#endif // FIXME
}
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index b2c1568a0cf6462dc771245ffc61aed51af184ec..c2f049277bac95b3ffda755b37d13bbeabad6539 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -45,13 +45,32 @@
template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
int _blockSize; // for block Matrix
- bool _isAllocated;
+ bool _isAllocated, _kspAllocated;
Mat _a;
Vec _b, _x;
+ KSP _ksp;
void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
+ void _kspCreate() {
+ _try(KSPCreate(PETSC_COMM_WORLD, &_ksp));
+ PC pc;
+ _try(KSPGetPC(_ksp, &pc));
+ // set some default options
+ _try(PCSetType(pc, PCILU));
+ _try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
+ _try(PCFactorSetLevels(pc, 1));
+ _try(KSPSetTolerances(_ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
+ // override the default options with the ones from the option
+ // database (if any)
+ _try(KSPSetFromOptions(_ksp));
+ _kspAllocated = true;
+ }
public:
- linearSystemPETSc(int blockSize = -1) : _isAllocated(false) {_blockSize = blockSize;}
- virtual ~linearSystemPETSc(){ clear(); }
+ linearSystemPETSc() : _isAllocated(false) {_blockSize = 0; _kspAllocated = false;}
+ virtual ~linearSystemPETSc(){
+ clear();
+ if (_kspAllocated)
+ _try (KSPDestroy (_ksp));
+ }
virtual bool isAllocated() const { return _isAllocated; }
virtual void allocate(int nbRows)
{
@@ -132,36 +151,39 @@ class linearSystemPETSc : public linearSystem<scalar> {
}
virtual void zeroMatrix()
{
- _try(MatZeroEntries(_a));
+ if (_isAllocated) {
+ _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
+ _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ _try(MatZeroEntries(_a));
+ }
}
virtual void zeroRightHandSide()
{
- _try(VecZeroEntries(_b));
+ if (_isAllocated) {
+ _try(VecAssemblyBegin(_b));
+ _try(VecAssemblyEnd(_b));
+ _try(VecZeroEntries(_b));
+ }
}
- virtual int systemSolve()
+ int systemSolve()
{
+ if (!_kspAllocated)
+ _kspCreate();
+ if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_sparsity")
+ _try(KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN));
+ else if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix")
+ _try(KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER));
+ else
+ _try(KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN));
_try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
_try(VecAssemblyBegin(_b));
_try(VecAssemblyEnd(_b));
- KSP ksp;
- _try(KSPCreate(PETSC_COMM_WORLD, &ksp));
- _try(KSPSetOperators(ksp, _a, _a, DIFFERENT_NONZERO_PATTERN));
- PC pc;
- _try(KSPGetPC(ksp, &pc));
- // set some default options
- _try(PCSetType(pc, PCILU));
- _try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
- _try(PCFactorSetLevels(pc, 1));
- _try(KSPSetTolerances(ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
- // override the default options with the ones from the option
- // database (if any)
- _try(KSPSetFromOptions(ksp));
- _try(KSPSolve(ksp, _b, _x));
- _try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
- PetscInt its;
- _try(KSPGetIterationNumber(ksp, &its));
- Msg::Info("%d iterations", its);
+ _try(KSPSolve(_ksp, _b, _x));
+ //_try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
+ //PetscInt its;
+ //_try(KSPGetIterationNumber(ksp, &its));
+ //Msg::Info("%d iterations", its);
return 1;
}
Mat &getMatrix(){ return _a; }