diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 0cbbb498d80e812e373a666cd72baea16d56c2d1..eccdc6bf93a18c5287aac8c55b7ab39bf017fe02 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1100,7 +1100,8 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value = myAssembler.getDofValue(v, 0, 1);
+ double value;
+ myAssembler.getDofValue(value, v, 0, 1);
std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if(itf == coordinates.end()){
SPoint3 p(0, 0, 0);
@@ -1348,8 +1349,9 @@ bool GFaceCompound::parametrize_conformal() const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value1 = myAssembler.getDofValue(v,0,1);
- double value2 = myAssembler.getDofValue(v,0,2);
+ double value1, value2;
+ myAssembler.getDofValue(value1, v,0,1);
+ myAssembler.getDofValue(value2, v,0,2);
coordinates[v] = SPoint3(value1,value2,0.0);
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 9deb77cb6bc77238fd811735a525e2a32737664c..7822b83d7ec5b8c598ebc8b526a24a7a3b76b538 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -371,8 +371,7 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
MVertex *v_3D = itv2->second;
- double value = myAssembler.getDofValue(v_3D, 0, 1);
- _sizes[v_2D] = value;
+ myAssembler.getDofValue(_sizes[v_2D], v_3D, 0, 1);
}
delete _lsys;
#endif
@@ -473,8 +472,7 @@ void backgroundMesh::updateSizes(GFace *_gf)
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
MVertex *v_3D = itv2->second;
- double value = myAssembler.getDofValue(v_3D, 0, 1);
- _sizes[v_2D] = value;
+ myAssembler.getDofValue(_sizes[v_2D], v_3D, 0, 1);
}
delete _lsys;
#endif
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 95957bee56a67395824e8432fd3c64457f43c155..281cde696c9e9b889d18c708d75074e97ffd82e5 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -599,8 +599,9 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
if ((*it)->onWhat()->dim() == 2){
SPoint2 param;
reparamMeshVertexOnFace((*it), gf, param);
- SPoint2 dparam (myAssembler.getDofValue((*it), 0, getTag()),
- myAssembler.getDofValue((*it), 1, getTag()));
+ SPoint2 dparam;
+ myAssembler.getDofValue(dparam[0], (*it), 0, getTag());
+ myAssembler.getDofValue(dparam[1], (*it), 1, getTag());
SPoint2 newp = param+dparam;
dx += newp.x() * newp.x() + newp.y() * newp.y();
(*it)->setParameter(0, newp.x());
@@ -715,9 +716,13 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
}
for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
- it->first->x() += myAssembler.getDofValue(it->first, 0, getTag());
- it->first->y() += myAssembler.getDofValue(it->first, 1, getTag());
- it->first->z() += myAssembler.getDofValue(it->first, 2, getTag());
+ double ax, ay, az;
+ myAssembler.getDofValue(ax, it->first, 0, getTag());
+ myAssembler.getDofValue(ay, it->first, 1, getTag());
+ myAssembler.getDofValue(az, it->first, 2, getTag());
+ it->first->x() += ax;
+ it->first->y() += ay;
+ it->first->z() += az;
}
// delete matrices and vectors
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 5d5a5f878e3c329cab904a3e63813c25a17f1b83..04630c7ae53d15b0cb8022530662ef3844a7934e 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -63,6 +63,8 @@ class fullVector
{
if(s == 0.)
for(int i = 0; i < _r; ++i) _data[i] = 0.;
+ else if (s == -1.)
+ for(int i = 0; i < _r; ++i) _data[i] = -_data[i];
else
for(int i = 0; i < _r; ++i) _data[i] *= s;
}
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index a5cb8df70e4480ece4de55d9fc82c069608dbb43..bb5634416319a46bdf538dbb9149f025e03760b3 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -188,7 +188,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
for(std::map<MVertex*,double >::iterator itv = distance_map2.begin();
itv !=distance_map2.end() ; ++itv){
MVertex *v = itv->first;
- double value = std::min(0.9999, myAssembler.getDofValue(v, 0, 1));
+ double value;
+ myAssembler.getDofValue(value, v, 0, 1);
+ value = std::min(0.9999, value);
double dist = -mu * log(1. - value);
itv->second = dist;
}
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index fbeea2c22e3a8dda43b855a75318918c3db6debf..d297e54425c6abb19e7586f38ef3622c6312d419 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -6,6 +6,7 @@
set(SRC
linearSystem.cpp
linearSystemCSR.cpp
+ linearSystemPETSc.cpp
groupOfElements.cpp
elasticityTerm.cpp
elasticitySolver.cpp
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 9df7b57f883e83a374654f98e8a15b5015f5d357..31a1fb782bc3b59d5b667a13840a382271b73a63 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -48,12 +48,18 @@ template<class T> struct dofTraits
{
typedef T VecType;
typedef T MatType;
+ inline static void mult (VecType &r, const MatType &m, const VecType &v) { r = m*v;}
+ inline static void neg (VecType &r) { r = -r;}
+ inline static void setToZero (VecType &r) { r = 0;}
};
-template<class T> struct dofTraits<fullVector<T> >
+template<class T> struct dofTraits<fullMatrix<T> >
{
- typedef fullVector<T> VecType;
+ typedef fullMatrix<T> VecType;
typedef fullMatrix<T> MatType;
+ inline static void mult (VecType &r, const MatType &m, const VecType &v) { m.mult(v, r);}
+ inline static void neg (VecType &r) { r.scale(-1.);}
+ inline static void setToZero (VecType &r) { r.scale(0);}
};
/*
template<> struct dofTraits<fullVector<std::complex<double> > >
@@ -159,24 +165,30 @@ class dofManager{
return dataVec(0.0);
}
- inline dataVec getDofValue(int ent, int type) const
+ inline void getDofValue(dataVec &v, 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;
+ if (it != fixed.end()){
+ v = it->second;
+ return;
+ }
}
{
std::map<Dof, int>::const_iterator it = unknown.find(key);
- if (it != unknown.end())
- return _current->getFromSolution(it->second);
+ if (it != unknown.end()){
+ v = _current->getFromSolution(it->second);
+ return;
+ }
}
- return dataVec(0.0);
+ dofTraits<T>::setToZero(v);
}
- inline dataVec getDofValue(MVertex *v, int iComp, int iField) const
+ inline dataVec getDofValue(dataVec &value, MVertex *v, int iComp, int iField) const
{
- return getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
+ getDofValue(value, 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());
@@ -189,8 +201,12 @@ class dofManager{
}
else{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
- if (itFixed != fixed.end()){
- _current->addToRightHandSide(itR->second, -value * itFixed->second);
+ if (itFixed != fixed.end()) {
+ // tmp = -value * itFixed->second
+ dataVec tmp(itFixed->second);
+ dofTraits<T>::mult(tmp , value, itFixed->second);
+ dofTraits<T>::neg(tmp);
+ _current->addToRightHandSide(itR->second, tmp);
}
}
}
@@ -220,7 +236,11 @@ class dofManager{
else{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C[j]);
if (itFixed != fixed.end()){
- _current->addToRightHandSide(NR[i], -m(i,j) * itFixed->second);
+ // tmp = -m(i,j) * itFixed->second
+ dataVec tmp(itFixed->second);
+ dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
+ dofTraits<T>::neg(tmp);
+ _current->addToRightHandSide(NR[i], tmp);
}
}
}
@@ -272,7 +292,11 @@ class dofManager{
{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(R[j]);
if (itFixed != fixed.end()){
- _current->addToRightHandSide(NR[i], -m(i, j) * itFixed->second);
+ // tmp = -m(i,j) * itFixed->second
+ dataVec tmp(itFixed->second);
+ dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
+ dofTraits<T>::neg(tmp);
+ _current->addToRightHandSide(NR[i], tmp);
}
}
}
@@ -336,5 +360,4 @@ class dofManager{
return 0;
}
};
-
#endif
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 047e83fdbc4a798b4e94e2deb8bb7acd471acc7f..65540b388c94fbea3910aa7c8b355ceeafd2d324 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -225,9 +225,9 @@ static double vonMises(dofManager<double> *a, MElement *e,
double valy[256];
double valz[256];
for (int k = 0; k < e->getNumVertices(); k++){
- valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
- valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
- valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
+ a->getDofValue(valx[k], e->getVertex(k), 0, _tag);
+ a->getDofValue(valy[k], e->getVertex(k), 1, _tag);
+ a->getDofValue(valz[k], e->getVertex(k), 2, _tag);
}
double gradux[3];
double graduy[3];
diff --git a/Solver/function.cpp b/Solver/function.cpp
index aadd3c0699dcaa41a3852188ed4074aadc17b99b..2d65bc040e931ddf15f1c3881bf1467e6348b0a7 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -447,7 +447,7 @@ void function::registerBindings(binding *b){
mb->setDescription("return the unique instance of this class");
cb->setParentClass<function>();
- cb = b->addClass<functionSolution>("functionNormals");
+ cb = b->addClass<functionNormals>("functionNormals");
cb->setDescription("A function to access the face normals. This is a single-instance class, use the 'get' member to access the instance.");
mb = cb->addMethod("get",&functionNormals::get);
mb->setDescription("return the unique instance of this class");
diff --git a/Solver/linearSystem.cpp b/Solver/linearSystem.cpp
index 32d0022012270307683079dfb7c76e92a6912e02..84193cd788e023e160c51d5d02d75012b06dd1cf 100644
--- a/Solver/linearSystem.cpp
+++ b/Solver/linearSystem.cpp
@@ -28,19 +28,13 @@ void linearSystem<double>::registerBindings(binding *b){
cb->setParentClass<linearSystem<double> >();
#endif
-#ifdef HAVE_PETSC
- cb = b->addClass<linearSystemPETSc<double> >("linearSystemPETSc");
- cb->setDescription("A linear system solver, based on PETSc");
- cm = cb->setConstructor<linearSystemPETSc<double> >();
- cm->setDescription ("A new PETSc<double> solver");
- cm->setArgNames(NULL);
- cb->setParentClass<linearSystem<double> >();
-#endif
-
cb = b->addClass<linearSystemFull<double> >("linearSystemFull");
cb->setDescription("A linear system solver, based on LAPACK (full matrices)");
cm = cb->setConstructor<linearSystemFull<double> >();
cm->setDescription ("A new Lapack based <double> solver");
cm->setArgNames(NULL);
cb->setParentClass<linearSystem<double> >();
+#ifdef HAVE_PETSC
+ linearSystemPETScRegisterBindings (b);
+#endif
}
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c77af8d464bb156544d4bc3cf5bde89ac1812911
--- /dev/null
+++ b/Solver/linearSystemPETSc.cpp
@@ -0,0 +1,110 @@
+#include "GmshConfig.h"
+#if defined HAVE_PETSC
+#include "linearSystemPETSc.h"
+#include "fullMatrix.h"
+template <>
+void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col, fullMatrix<double> val)
+{
+ for (int ii = 0; ii<val.size1(); ii++)
+ for (int jj = 0; jj < ii; jj ++) {
+ double buff = val(ii,jj);
+ val(ii,jj) = val (jj,ii);
+ val(jj,ii) = buff;
+ }
+ PetscInt i = row, j = col;
+ _try(MatSetValuesBlocked(_a, 1, &i, 1, &j, &val(0,0), ADD_VALUES));
+}
+template<>
+void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row, fullMatrix<double> val)
+{
+ for (int ii = 0; ii < _blockSize; ii ++) {
+ PetscInt i = row * _blockSize + ii;
+ _try(VecSetValues(_b, 1, &i, &val(ii,0), ADD_VALUES));
+ }
+}
+
+template<>
+fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromMatrix(int row, int col) const
+{
+ Msg::Error("getFromMatrix not implemented for PETSc");
+ return fullMatrix<double>(0,0);
+}
+
+template<>
+fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide(int row) const
+{
+ fullMatrix<double> val(_blockSize,1);
+ PetscScalar *tmp;
+ _try(VecGetArray(_b, &tmp));
+ for (int i=0; i<_blockSize; i++) {
+ PetscScalar s = tmp[row*_blockSize+i];
+#if defined(PETSC_USE_COMPLEX)
+ val(i,0) = s.real();
+#else
+ val(i,0) = s;
+#endif
+ }
+ _try(VecRestoreArray(_b, &tmp));
+ return val;
+}
+
+template<>
+fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row) const
+{
+ fullMatrix<double> val(_blockSize,1);
+ PetscScalar *tmp;
+ _try(VecGetArray(_x, &tmp));
+ for (int i=0; i<_blockSize; i++) {
+ PetscScalar s = tmp[row*_blockSize+i];
+#if defined(PETSC_USE_COMPLEX)
+ val(i,0) = s.real();
+#else
+ val(i,0) = s;
+#endif
+ }
+ _try(VecRestoreArray(_x, &tmp));
+ return val;
+}
+
+template<>
+void linearSystemPETSc<fullMatrix<double> >::allocate(int nbRows) {
+ clear();
+ _try(MatCreate(PETSC_COMM_WORLD, &_a));
+ _try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows*_blockSize, nbRows*_blockSize));
+ _try(MatSetType(_a, MATSEQBAIJ));
+ // override the default options with the ones from the option
+ // database (if any)
+ _try(MatSetFromOptions(_a));
+ _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 4, NULL)); //todo preAllocate off-diagonal part
+ //_try(MatMPIBAIJSetPreallocation(_a, _blockSize, 4, NULL, 0, NULL)); //todo preAllocate off-diagonal part
+ _try(VecCreate(PETSC_COMM_WORLD, &_x));
+ _try(VecSetSizes(_x, PETSC_DECIDE, nbRows*_blockSize));
+ // override the default options with the ones from the option
+ // database (if any)
+ _try(VecSetFromOptions(_x));
+ _try(VecDuplicate(_x, &_b));
+ _isAllocated = true;
+}
+
+#include "Bindings.h"
+void linearSystemPETScRegisterBindings(binding *b) {
+ classBinding *cb;
+ methodBinding *cm;
+ cb = b->addClass<linearSystemPETSc<double> >("linearSystemPETSc");
+ cb->setDescription("A linear system solver, based on PETSc");
+ cm = cb->setConstructor<linearSystemPETSc<double> >();
+ cm->setDescription ("A new PETSc<double> solver");
+ cb->setParentClass<linearSystem<double> >();
+ cm->setArgNames(NULL);
+ cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<double> >::systemSolve);
+ cm->setDescription("compute x = A^{-1}b");
+
+ cb = b->addClass<linearSystemPETSc<fullMatrix<double> > >("linearSystemPETScBlock");
+ cb->setDescription("A linear system solver, based on PETSc");
+ cm = cb->setConstructor<linearSystemPETSc<fullMatrix<double> >, int>();
+ cm->setDescription ("A new PETScBlock<double> solver (we probably should get rid of the blockSize argument)");
+ cm->setArgNames("blockSize", NULL);
+ cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<double> >::systemSolve);
+ cm->setDescription("compute x = A^{-1}b");
+}
+#endif // HAVE_PETSC
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index acc4a040e0147c60b9f0c25dbd52ff936c36508c..d17a14c333e13d8dd001a04007fa240ac19f6806 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -44,13 +44,13 @@
template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
- private:
+ int _blockSize; // for block Matrix
bool _isAllocated;
Mat _a;
Vec _b, _x;
void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
public:
- linearSystemPETSc() : _isAllocated(false) {}
+ linearSystemPETSc(int blockSize = -1) : _isAllocated(false) {_blockSize = blockSize;}
virtual ~linearSystemPETSc(){ clear(); }
virtual bool isAllocated() const { return _isAllocated; }
virtual void allocate(int nbRows)
@@ -83,12 +83,6 @@ class linearSystemPETSc : public linearSystem<scalar> {
}
_isAllocated = false;
}
- virtual void addToMatrix(int row, int col, scalar val)
- {
- PetscInt i = row, j = col;
- PetscScalar s = val;
- _try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
- }
virtual scalar getFromMatrix(int row, int col) const
{
Msg::Error("getFromMatrix not implemented for PETSc");
@@ -113,6 +107,12 @@ class linearSystemPETSc : public linearSystem<scalar> {
return s;
#endif
}
+ virtual void addToMatrix(int row, int col, scalar val)
+ {
+ PetscInt i = row, j = col;
+ PetscScalar s = val;
+ _try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
+ }
virtual scalar getFromSolution(int row) const
{
PetscScalar *tmp;
@@ -152,6 +152,9 @@ class linearSystemPETSc : public linearSystem<scalar> {
// override the default options with the ones from the option
// database (if any)
_try(KSPSetFromOptions(ksp));
+ PetscViewer view;
+ PetscViewerASCIIOpen(PETSC_COMM_WORLD,"mat.out", &view);
+ _try(MatView(_a, view));
_try(KSPSolve(ksp, _b, _x));
_try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
PetscInt its;
@@ -161,6 +164,10 @@ class linearSystemPETSc : public linearSystem<scalar> {
Mat &getMatrix(){ return _a; }
};
+class binding;
+void linearSystemPETScRegisterBindings(binding *b);
+
+
#else
template <class scalar>
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 73637c6874988773357f77d50ec26b078ab3591a..29acbb5a29e68499d37a2b2866448167afabd3b0 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -1199,7 +1199,8 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
int count = 0;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value = myAssembler.getDofValue(v, 0, 1);
+ double value;
+ myAssembler.getDofValue(value, v, 0, 1);
if (step == 0)solution[v] = SPoint2(value,0);
else solution[v] = SPoint2(solution[v][0],value);
}