diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index 4707da41035bc78fd344d077815f64deede0d0dc..a1470c97c72adb8d6ee4547eae23eef39e7108a0 100644 --- a/Geo/GFaceCompound.cpp +++ b/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; diff --git a/Solver/diagBCTerm.h b/Solver/diagBCTerm.h index ce59034a8b006070e572288f9f499db9e50ca2e4..6e13dfd2f30d5b828fc435f3bc3f1f2d5d25a5e0 100644 --- a/Solver/diagBCTerm.h +++ b/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; } } diff --git a/Solver/dofManager.h b/Solver/dofManager.h index 5e0f0605728a6cd05eba0a7b6698aba175f3f6fa..41f39ef10479d75f43a92f6d64d3657d9c601ce6 100644 --- a/Solver/dofManager.h +++ b/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) { diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp index f162a08f5c1b785bd0dbfe77da400d6aff97eca1..d87b8d97621f68aa77d2783d6c24ee6a4ab363aa 100644 --- a/Solver/eigenSolver.cpp +++ b/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));