verbosity=0; load "iovtk" load "gmsh" load "msh3" func int Mylab(int i){ return i; } mesh TH=square(2,2); plot(TH, wait=1); int nbtri=TH.nt; //number of coarse elements in the coarse mesh fespace Tri(TH,P0); // P0 on coarse mesh Tri ChiK; for (int k=0; k0.1, split=1); // a mesh made of only one triangle; ChiK[][k]=0; // P0 function used to mark the element i int m=0; THK = change(THK, flabel=Mylab(m++)); mesh ThK = trunc(THK, 1, split=100); fespace VhK(ThK, P1); //FE space to solve the local problem VhK[int] PhiK(4); //FE Array to temporary storage of the 4 basis functions per coarse element (1 bubble + 3 edges) fespace VHMsFEMK(THK,[P0,P0edge]); // FE MsFEM on coarse mesh to compute velocity varf VZEROK(u,v)=int2d(THK)(0.*u*v); matrix STIFFK = VZEROK(VHMsFEMK,VHMsFEMK); matrix MASSK = VZEROK(VHMsFEMK,VHMsFEMK); string basisstorePhi ="./BASIS_TEST/Basis_ELEM_"+k+".txt"; { ifstream readbasisPhi(basisstorePhi); for (int j=0; j<4; j++){ readbasisPhi >> PhiK[j][]; } //loading of the 4 MsFEM basis functions per edge } //******************Build Variational formulations for local stiffness matrix*************************************************** varf a(u,v)= int2d(ThK, qforder=6)(dx(u)*dx(v)+dy(u)*dy(v)); varf b(u,v)= int2d(ThK, qforder=6)(u*v); //******************Build the matrix for local stiffness matrix*************************************************** matrix AK = a(VhK,VhK); // Matrix to compute componument S_{ij} of stiffness matrices matrix BK = b(VhK, VhK); // Matrix to compute componument M_{ij} of mass matrices for (int j=0; j<4; j++) { int I=VHMsFEMK(0,j); for (int l=0; l<4; l++){ int J=VHMsFEMK(0,l); real[int] AKinter=AK*PhiK[l][]; real AKMsFEM= PhiK[j][]'*AKinter; STIFFK(I,J)=STIFFK(I,J)+AKMsFEM;}} for (int j=0; j<4; j++) { int I=VHMsFEMK(0,j); for (int l=0; l<4; l++){ int J=VHMsFEMK(0,l); real[int] BKinter=BK*PhiK[l][];; real BKMsFEM=PhiK[j][]'*BKinter; MASSK(I,J)=MASSK(I,J)+BKMsFEM;}} int nev=2; real sigma=0; matrix OP = MASSK - sigma * STIFFK; real[int] ev(nev); //to store the nev eigenvalue VHMsFEMK[int] eV(nev); //to store the nev eigenvector l = EigenValue(OP, STIFFK, sym=true, sigma=sigma, value=ev, vector=eV,tol=1e-10, maxit=0, ncv=0); for (int i = 0; i < l; i++){ cout << "lambda[" << i << "] = " << ev[i] << endl;} }