From 1b002f74d61523751b2bc2b9cf93905a0b05d602 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 19 Feb 2010 10:06:24 +0000
Subject: [PATCH]
---
Geo/GFaceCompound.cpp | 15 ++++++++++-----
Solver/diagBCTerm.h | 7 +++----
Solver/dofManager.h | 2 --
Solver/eigenSolver.cpp | 24 ++++++++++++++++++++++--
4 files changed, 35 insertions(+), 13 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 4707da4103..a1470c97c7 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 ce59034a8b..6e13dfd2f3 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 5e0f060572..41f39ef104 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 f162a08f5c..d87b8d9762 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));
--
GitLab