diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a1470c97c72adb8d6ee4547eae23eef39e7108a0..ae5756a6b50ac6584c4a1f9a8928646e967abcb8 100644
--- a/Geo/GFaceCompound.cpp
+++ b/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,48 +1138,91 @@ 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());
-
- int k = 0;
- std::vector<std::complex<double> > &ev = eig.getEigenVector(0);
- 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;
+ 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){
+ MVertex *v = *itv;
+ double paramu = ev[k].real();
+ double paramv = ev[k+1].real();
+ coordinates[v] = SPoint3(paramu,paramv,0.0);
+ 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();
+
+ return checkFolding(ordered);
+
}
-
- lsysA->clear();
- lsysB->clear();
-
- //check for folding
- return checkFolding(ordered);
+ else return false;
#else
return false;
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index d87b8d97621f68aa77d2783d6c24ee6a4ab363aa..9a582c7e564f905a55ad4c4292498b203481ae4c 100644
--- a/Solver/eigenSolver.cpp
+++ b/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,15 +87,18 @@ void eigenSolver::solve(int numEigenValues, std::string which)
// solve
Msg::Info("SLEPc solving...");
+ double t1 = Cpu();
_try(EPSSolve(eps));
-
+
// check convergence
int its;
_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)
@@ -110,55 +107,70 @@ void eigenSolver::solve(int numEigenValues, std::string which)
Msg::Error("The operator is nonsymmetric");
// get number of converged approximate eigenpairs
- PetscInt nconv;
- _try(EPSGetConverged(eps, &nconv));
- Msg::Info("SLEPc number of converged eigenpairs: %d", nconv);
-
- // ignore additional eigenvalues if we get more than what we asked
- if(nconv > nev) nconv = nev;
-
- Vec xr, xi;
- _try(MatGetVecs(A, PETSC_NULL, &xr));
- _try(MatGetVecs(A, PETSC_NULL, &xi));
- Msg::Info(" Re[EigenValue] Im[EigenValue]"
- " Relative error");
- for (int i = 0; i < nconv; i++){
- PetscScalar kr, ki;
- _try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi));
- PetscReal error;
- _try(EPSComputeRelativeError(eps, i, &error));
+ PetscInt nconv;
+ _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;
+
+ Vec xr, xi;
+ _try(MatGetVecs(A, PETSC_NULL, &xr));
+ _try(MatGetVecs(A, PETSC_NULL, &xi));
+ Msg::Info(" Re[EigenValue] Im[EigenValue]"
+ " Relative error");
+ for (int i = 0; i < nconv; i++){
+ PetscScalar kr, ki;
+ _try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi));
+ PetscReal error;
+ _try(EPSComputeRelativeError(eps, i, &error));
#if defined(PETSC_USE_COMPLEX)
- PetscReal re = PetscRealPart(kr);
- PetscReal im = PetscImaginaryPart(kr);
+ PetscReal re = PetscRealPart(kr);
+ PetscReal im = PetscImaginaryPart(kr);
#else
- PetscReal re = kr;
- PetscReal im = ki;
+ PetscReal re = kr;
+ PetscReal im = ki;
#endif
- Msg::Info("EIG %03d %s%.16e %s%.16e %3.6e",
- i, (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
-
- // store eigenvalues and eigenvectors
- _eigenValues.push_back(std::complex<double>(re, im));
- PetscScalar *tmpr, *tmpi;
- _try(VecGetArray(xr, &tmpr));
- _try(VecGetArray(xi, &tmpi));
- std::vector<std::complex<double> > ev(N);
- for(int i = 0; i < N; i++){
+ Msg::Info("EIG %03d %s%.16e %s%.16e %3.6e",
+ i, (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
+
+ // store eigenvalues and eigenvectors
+ _eigenValues.push_back(std::complex<double>(re, im));
+ PetscScalar *tmpr, *tmpi;
+ _try(VecGetArray(xr, &tmpr));
+ _try(VecGetArray(xi, &tmpi));
+ std::vector<std::complex<double> > ev(N);
+ for(int i = 0; i < N; i++){
#if defined(PETSC_USE_COMPLEX)
- ev[i] = tmpr[i];
+ ev[i] = tmpr[i];
#else
- ev[i] = std::complex<double>(tmpr[i], tmpi[i]);
+ ev[i] = std::complex<double>(tmpr[i], tmpi[i]);
#endif
+ }
+ _eigenVectors.push_back(ev);
}
- _eigenVectors.push_back(ev);
+
+ // cleanup
+ _try(EPSDestroy(eps));
+ _try(VecDestroy(xr));
+ _try(VecDestroy(xi));
+ _try(SlepcFinalize());
+
}
+
+
- // cleanup
- _try(EPSDestroy(eps));
- _try(VecDestroy(xr));
- _try(VecDestroy(xi));
- _try(SlepcFinalize());
- Msg::Info("SLEPc done");
+ if(reason == EPS_CONVERGED_TOL){
+ Msg::Info("SLEPc done");
+ return true;
+ }
+ else{
+ Msg::Info("SLEPc failed");
+ return false;
+ }
+
}
#endif
diff --git a/Solver/eigenSolver.h b/Solver/eigenSolver.h
index a0ddcd0b9b8fb46af2a477a7729d64fc38475bf9..d4c9394d81111199c364b12c2d0db27cc28dec4f 100644
--- a/Solver/eigenSolver.h
+++ b/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");