block linear solver for dg

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);
   }
 

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

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

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;
   }

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;
   }

Solver/CMakeLists.txt

@@ -6,6 +6,7 @@
 set(SRC
   linearSystem.cpp
   linearSystemCSR.cpp
+  linearSystemPETSc.cpp
   groupOfElements.cpp
   elasticityTerm.cpp
   elasticitySolver.cpp

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);
     }
-  inline dataVec getDofValue(MVertex *v, int iComp, int iField) const
+    dofTraits<T>::setToZero(v);
+  }
+  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());
@@ -190,7 +202,11 @@ class dofManager{
       else{
         typename std::map<Dof,  dataVec>::iterator itFixed = fixed.find(C);
         if (itFixed != fixed.end()) {
-          _current->addToRightHandSide(itR->second, -value * itFixed->second);
+          // 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

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];

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");

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
 }

Solver/linearSystemPETSc.cpp

+#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

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>

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);
     }