Temporal stokes ok

Numeric/fullMatrix.h

@@ -123,6 +123,14 @@ class fullMatrix
     #endif
     (*this)(r, c) = v; 
   }
+  inline scalar norm()
+  {
+    scalar n = 0.;
+    for(int i = 0; i < _r; ++i) 
+      for(int j = 0; j < _c; ++j) 
+	n += (*this)(i, j) * (*this)(i, j);
+    return sqrt(n);
+  }
   fullMatrix(scalar *original, int r, int c)
   {
     _r = r;

Solver/linearSystemPETSc.cpp

@@ -82,6 +82,8 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
   if (_blockSize == 0)
     Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
   clear();
+  //printf("allocate linear system petsc \n");
+  //_try(PetscOptionsInsertString("-ksp_monitor_true_residual -ksp_rtol 1e-10"));
   _try(MatCreate(PETSC_COMM_WORLD, &_a)); 
   _try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows * _blockSize, nbRows * _blockSize));
   _try(MatSetType(_a, MATSEQBAIJ));

Solver/linearSystemPETSc.h

@@ -55,7 +55,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
     PC pc;
     _try(KSPGetPC(_ksp, &pc));
     // set some default options
-    _try(PCSetType(pc, PCILU));
+    _try(PCSetType(pc, PCLU));//LU for direct solver and PCILU for indirect solver
     _try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
     _try(PCFactorSetLevels(pc, 1));
     _try(KSPSetTolerances(_ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));