No commit message

No commit message

Geo/GFaceCompound.cpp

@@ -1124,6 +1124,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
     MVertex *v = *itv;
     myAssembler.numberVertex(v, 0, 1);
     myAssembler.numberVertex(v, 0, 2);
+
   }
 
   simpleFunction<double> ONE(1.0);
@@ -1151,18 +1152,22 @@ bool GFaceCompound::parametrize_conformal_spectral() const
      myAssembler.numberVertex(v, 0, 1);
      myAssembler.numberVertex(v, 0, 2);
    }
-   diagBCTerm diag(0, 1, &ONE);
+   
+   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]);
-       diag.addToMatrix(myAssembler, &se);
+       diag1.addToMatrix(myAssembler, &se);
+       diag2.addToMatrix(myAssembler, &se);
      }
    }
-   myAssembler.setCurrentMatrix("A");
-  //-------------------------------
+
+   //-------------------------------
    eigenSolver eig(&myAssembler, "A" ); //, "B");
-   eig.solve(2, "smallest");
+   //eig.solve(1, "largest");
+   eig.solve(1, "smallestReal");
    //printf("num eigenvalues =%d \n", eig.getNumEigenValues());
    
    int k = 0;

Solver/diagBCTerm.h

@@ -48,10 +48,9 @@ class diagBCTerm : public femTerm<double> {
 	m(j,k) = 0.0;
 	m(k,j) = 0.0;
       }
-      m(j,j)  = 1.0;
-/*       MVertex *v = e->getVertex(j); */
-/*       if( v->onWhat()->dim() < 2 )  m(j,j) = (nbNodes - 1) ; */
-/*       else m(j,j) = 0.0; */
+      MVertex *v = e->getVertex(j);
+      if( v->onWhat()->dim() < 2 ) m(j,j) = 1.0; 
+      else m(j,j) = 0.0;
     }
   } 
   

Solver/dofManager.h

@@ -105,8 +105,6 @@ class dofManager{
  dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2) : _current(l1) { 
     _linearSystems.insert(std::make_pair("A", l1));
     _linearSystems.insert(std::make_pair("B", l2));
-    //_linearSystems.["A"] = l1; 
-    //_linearSystems["B"] = l2;
   }
   inline void fixDof(Dof key, const dataVec &value)
   {

Solver/eigenSolver.cpp

@@ -32,7 +32,14 @@ void eigenSolver::solve(int numEigenValues, std::string which)
   _try(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
   PetscInt N, M;
   _try(MatGetSize(A, &N, &M));
- 
+
+  PetscInt N2, M2;
+  if (_B) {
+    _try(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
+    _try(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
+    _try(MatGetSize(B, &N2, &M2));
+  } 
+
   // generalized eigenvalue problem A x - \lambda B x = 0
   EPS eps;
   _try(EPSCreate(PETSC_COMM_WORLD, &eps));
@@ -43,6 +50,18 @@ void eigenSolver::solve(int numEigenValues, std::string which)
   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));
 
@@ -54,10 +73,11 @@ void eigenSolver::solve(int numEigenValues, std::string which)
     _try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
   if(which == "smallest")
     _try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE));
+  else if(which == "smallestReal")
+    _try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL));
   else if(which == "largest")
     _try(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE));
 
-
   // print info
   const EPSType type;
   _try(EPSGetType(eps, &type));