spectral parametrization with slepc

Geo/GFaceCompound.cpp

@@ -505,6 +505,7 @@ bool GFaceCompound::parametrize() const
     fillNeumannBCS();
     bool withoutFolding = parametrize_conformal_spectral() ;
     printStuff();
+    //exit(1);
     if ( withoutFolding == false ){
       Msg::Warning("$$$ Parametrization switched to harmonic map");
       parametrize(ITERU,HARMONIC); 
@@ -530,10 +531,6 @@ bool GFaceCompound::parametrize() const
     buildOct();
   }
 
-
-
-
-
   if (checkAspectRatio() > AR_MAX){
     Msg::Warning("Geometrical aspect ratio too high");
     //exit(1);
@@ -1119,12 +1116,10 @@ bool GFaceCompound::parametrize_conformal_spectral() const
 
   //-------------------------------
   myAssembler.setCurrentMatrix("A");
-
   for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
     MVertex *v = *itv;
     myAssembler.numberVertex(v, 0, 1);
     myAssembler.numberVertex(v, 0, 2);
-
   }
 
   simpleFunction<double> ONE(1.0);
@@ -1143,33 +1138,59 @@ bool GFaceCompound::parametrize_conformal_spectral() const
       cross21.addToMatrix(myAssembler, &se);
     }
   }
+  double epsilon = 1.e-6;
+  for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+      MVertex *v = *itv;
+      if (std::find(ordered.begin(), ordered.end(), v) == ordered.end() ){
+        myAssembler.assemble(v, 0, 1, v, 0, 1,  epsilon);
+        myAssembler.assemble(v, 0, 2, v, 0, 2,  epsilon);
+      }
+  }
 
   //-------------------------------
    myAssembler.setCurrentMatrix("B");
-
    for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
      MVertex *v = *itv;
      myAssembler.numberVertex(v, 0, 1);
      myAssembler.numberVertex(v, 0, 2);
    }
       
-   diagBCTerm diag1(0, 1, &ONE);
-   diagBCTerm diag2(0, 2, &ONE);
-   it = _compound.begin(); 
-   for( ; it != _compound.end() ; ++it){
-     for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
-       SElement se((*it)->triangles[i]);
-       diag1.addToMatrix(myAssembler, &se);
-       diag2.addToMatrix(myAssembler, &se);
+   double small = 0.0;
+   for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+     MVertex *v = *itv;
+     if (std::find(ordered.begin(), ordered.end(), v) == ordered.end() ){
+       myAssembler.assemble(v, 0, 1, v, 0, 1,  small);
+       myAssembler.assemble(v, 0, 2, v, 0, 2,  small);
+     }
+     else{
+       myAssembler.assemble(v, 0, 1, v, 0, 1,  1.0);
+       myAssembler.assemble(v, 0, 2, v, 0, 2,  1.0);
      }
    }
+//    int NB = ordered.size();
+//    for(std::vector<MVertex *>::iterator itv1 = ordered.begin(); itv1 !=ordered.end() ; ++itv1){
+//      for(std::vector<MVertex *>::iterator itv2 = ordered.begin(); itv2 !=ordered.end() ; ++itv2){
+//        myAssembler.assemble(*itv1, 0, 1, *itv2, 0, 1,  -1/NB);
+//        myAssembler.assemble(*itv1, 0, 2, *itv2, 0, 2,  -1/NB);
+//      }
+//    }
+
+//    diagBCTerm diag1(0, 1, &ONE);
+//    diagBCTerm diag2(0, 2, &ONE);
+//    it = _compound.begin(); 
+//    for( ; it != _compound.end() ; ++it){
+//      for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+//        SElement se((*it)->triangles[i]);
+//        diag1.addToMatrix(myAssembler, &se);
+//        diag2.addToMatrix(myAssembler, &se);
+//      }
+//    }
 
    //-------------------------------
-   eigenSolver eig(&myAssembler, "A" ); //, "B");
-   //eig.solve(1, "largest");
-   eig.solve(1, "smallestReal");
-   //printf("num eigenvalues =%d \n", eig.getNumEigenValues());
+   eigenSolver eig(&myAssembler, "B" , "A", true);
+   bool converged = eig.solve(2, "largest");
      
+   if(converged) {
      int k = 0;
      std::vector<std::complex<double> > &ev = eig.getEigenVector(0); 
      for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -1180,12 +1201,29 @@ bool GFaceCompound::parametrize_conformal_spectral() const
        k = k+2;
      }
      
+     //if folding take second sallest eigenvalue
+     bool noFolding = checkFolding(ordered);
+     if (!noFolding ){
+       coordinates.clear();
+       int k = 0;
+       std::vector<std::complex<double> > &ev = eig.getEigenVector(1); 
+       for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+	 MVertex *v = *itv;
+	 double paramu = ev[k].real();
+	 double paramv = ev[k+1].real();
+	 coordinates[v] = SPoint3(paramu,paramv,0.0);
+	 k = k+2;
+       }
+     }
+
      lsysA->clear();
      lsysB->clear();
 
-   //check for folding
      return checkFolding(ordered);
 
+   }
+   else return false;
+
 #else
    return false;
 

Solver/eigenSolver.cpp

@@ -4,13 +4,14 @@
 // bugs and problems to <gmsh@geuz.org>.
 
 #include "eigenSolver.h"
+#include "OS.h"
 
 #if defined(HAVE_SLEPC)
 
 #include <slepceps.h>
 
 eigenSolver::eigenSolver(dofManager<double> *manager, std::string A, 
-                         std::string B) : _A(0), _B(0)
+                         std::string B, bool hermitian) : _A(0), _B(0),_hermitian(hermitian)
 {
   if(A.size()){
     _A = dynamic_cast<linearSystemPETSc<double>*>(manager->getLinearSystem(A));
@@ -20,11 +21,13 @@ eigenSolver::eigenSolver(dofManager<double> *manager, std::string A,
     _B = dynamic_cast<linearSystemPETSc<double>*>(manager->getLinearSystem(B));
     if(!_B) Msg::Error("Could not find PETSc system '%s'", B.c_str());
   }
+
 }
 
-void eigenSolver::solve(int numEigenValues, std::string which)
+bool eigenSolver::solve(int numEigenValues, std::string which)
 {
-  if(!_A) return;
+
+  if(!_A) return false;
   Mat A = _A->getMatrix();
   Mat B = _B ? _B->getMatrix() : PETSC_NULL;
    
@@ -44,26 +47,17 @@ void eigenSolver::solve(int numEigenValues, std::string which)
   EPS eps;
   _try(EPSCreate(PETSC_COMM_WORLD, &eps));
   _try(EPSSetOperators(eps, A, B));
-  bool hermitian = false; // FIXME
-  if(hermitian)
+  if(_hermitian)
     _try(EPSSetProblemType(eps, _B ? EPS_GHEP : EPS_HEP));
   else
     _try(EPSSetProblemType(eps, _B ? EPS_GNHEP : EPS_NHEP));
 
-//   EPSSetProblemType(eps, EPS_GNHEP);
-//   ST st;
-//   EPSGetST(eps,&st);
-//   KSP ksp;
-//   PC pc;
-//   STGetKSP(st,&ksp); 
-//   KSPGetPC(ksp, &pc);
-//   PCSetType(pc, PCJACOBI);
-//   PCFactorSetMatOrderingType(pc, MATORDERING_RCM);
-//   STSetType(st,STSINV);
-//   STSetShift(st,0.0);
-
   // set some default options
-  _try(EPSSetDimensions(eps, 5, PETSC_DECIDE, PETSC_DECIDE));
+  _try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
+  _try(EPSSetTolerances(eps,1.e-7,100));
+  _try(EPSSetType(eps,EPSARNOLDI)); //EPSKRYLOVSCHUR is default
+  //_try(EPSSetType(eps,EPSARPACK)); 
+  //_try(EPSSetType(eps,EPSPOWER));
 
   // override these options at runtime, petsc-style
   _try(EPSSetFromOptions(eps));
@@ -93,6 +87,7 @@ void eigenSolver::solve(int numEigenValues, std::string which)
 
   // solve
   Msg::Info("SLEPc solving...");
+  double t1 = Cpu();
   _try(EPSSolve(eps));
  
   // check convergence
@@ -100,8 +95,10 @@ void eigenSolver::solve(int numEigenValues, std::string which)
   _try(EPSGetIterationNumber(eps, &its));
   EPSConvergedReason reason;
   _try(EPSGetConvergedReason(eps, &reason));
-  if(reason == EPS_CONVERGED_TOL)
-    Msg::Info("SLEPc converged in %d iterations", its);
+  if(reason == EPS_CONVERGED_TOL){
+    double t2 = Cpu();
+    Msg::Info("SLEPc converged in %d iterations (%g s)", its, t2-t1);
+  }
   else if(reason == EPS_DIVERGED_ITS)
     Msg::Error("SLEPc diverged after %d iterations", its);
   else if(reason == EPS_DIVERGED_BREAKDOWN)
@@ -114,6 +111,8 @@ void eigenSolver::solve(int numEigenValues, std::string which)
     _try(EPSGetConverged(eps, &nconv));
     Msg::Info("SLEPc number of converged eigenpairs: %d", nconv);
     
+    if (nconv>0) {
+
     // ignore additional eigenvalues if we get more than what we asked
     if(nconv > nev) nconv = nev;
     
@@ -158,7 +157,20 @@ void eigenSolver::solve(int numEigenValues, std::string which)
     _try(VecDestroy(xr));
     _try(VecDestroy(xi));
     _try(SlepcFinalize());
+
+  }
+
+
+  
+  if(reason == EPS_CONVERGED_TOL){
     Msg::Info("SLEPc done");
+    return true;
+  }
+  else{
+    Msg::Info("SLEPc failed");
+    return false;
+  }
+
 }
 
 #endif

Solver/eigenSolver.h

@@ -20,13 +20,14 @@
 class eigenSolver{
  private:
   linearSystemPETSc<double> *_A, *_B;
+  bool _hermitian;
   std::vector<std::complex<double> > _eigenValues;
   std::vector<std::vector<std::complex<double> > > _eigenVectors;
   void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
  public:
   eigenSolver(dofManager<double> *manager, std::string A, 
-              std::string B="");
-  void solve(int numEigenValues=0, std::string which="");
+              std::string B="", bool hermitian=false);
+  bool solve(int numEigenValues=0, std::string which="");
   int getNumEigenValues(){ return _eigenValues.size(); }
   std::complex<double> getEigenValue(int num){ return _eigenValues[num]; }
   std::vector<std::complex<double> > &getEigenVector(int num){ return _eigenVectors[num]; }
@@ -39,7 +40,7 @@ class eigenSolver{
   std::vector<std::complex<double> > _dummy;
  public:
   eigenSolver(dofManager<double> *manager, std::string A, 
-              std::string B=""){}
+              std::string B="", bool hermitian=false){}
   void solve(int numEigenValues=0, std::string which="")
   {
     Msg::Error("Eigen solver requires SLEPc");