#include "xlife++.h" #include #include using namespace xlifepp; using namespace std; int main(int argc, char** argv) { init(_lang=en); //mandatory initialization of xlifepp //Change the value of choice to choose between uniform and Gauss-Lobatto interpolation //The result will be in a file name degCond_UI or deg_Cond_GL //In this file on each ligne you have degree CN_S CN_M //Where CN_S is the condition number of the stiffness matrix and CN_M is the condition number of the mass matrix //choice(1) for uniform interpolation //choice(0) for Gauss-Lobatto interpolation int choice(1); Int nbvp(6); //number of eigenvalue computed Int dMin(4); //degree min Int dMax(20); //degree max String name; FESubType type; if(choice){ name = "UI"; type = _standard; } else{ name = "GL"; type = _GaussLobatto; } //Mesh and domain Square sq(_origin=Point(0.,0.),_length=pi_,_nnodes=2,_domain_name="Omega",_side_names="Gamma"); Mesh mesh(sq,_quadrangle,1,_structured); Domain Omega = mesh.domain("Omega"); Domain Gamma = mesh.domain("Gamma"); std::ofstream myfile; myfile.precision(17); String fileName("degCond_"+name); myfile.open(fileName.c_str()); for(int d=dMin; d<=dMax; d++){ cout << "---------------------------" << endl; cout << "degree = " << d << endl; //Space and unknown Space Vd(Omega,interpolation(Lagrange,type,d,H1),"Vd"); Unknown u(Vd, "u"); TestFunction v(u, "v"); //Bilinear forme BilinearForm auv = intg(Omega,grad(u)|grad(v)); BilinearForm buv = intg(Omega,u*v); EssentialConditions ecs = ((u|Gamma) = 0); //TermMatrix and Solver cpuTime(); TermMatrix A(auv, ecs, ReductionMethod(_pseudoReduction,10.), "A"); TermMatrix B(buv, ecs, ReductionMethod(_pseudoReduction,0.1), "B"); cpuTime("Time TermMatrix"); EigenElements eAmax = eigenSolve(A,_nev=nbvp,_which="LM"); EigenElements eAmin = eigenSolve(A,_nev=nbvp,_sigma=-1.); EigenElements eBmax = eigenSolve(B,_nev=nbvp,_which="LM"); EigenElements eBmin = eigenSolve(B,_nev=nbvp,_sigma=0.); cpuTime("Time eigenSolve"); myfile << d << " "; myfile << eAmax.value(nbvp).real()/eAmin.value(1).real() << " "; myfile << eBmax.value(nbvp).real()/eBmin.value(1).real() << endl; } myfile.close(); cout << "#### Fini ####" << endl; return 0; }