1 /*******************************************************
 2  *  This file is part of the FreeFEM tutorial made by  *
 3  *      Pierre Jolivet <pierre@joliv.et>               *
 4  *                                                     *
 5  *  See https://joliv.et/FreeFem-tutorial for more     *
 6  *                                                     *
 7  *   Description: eigenvectors of the                  *
 8  *                Steklov–Poincaré operator            *
 9  *******************************************************/
 
11 load "PETSc"
12 macro dimension()3// EOM
13 include "macro_ddm.idp"
 
15 include "cube.idp"
16 macro grad(u)[dx(u), dy(u), dz(u)]// EOM
 
18 int[int] LL = [1,2, 1,1, 1,1];
19 mesh3 cb = cube(20, 20, 20, [x, y, z], label = LL);
20 Mat A;
21 func Pk = P2;
22 MatCreate(cb, A, Pk);
23 fespace Vh(cb, Pk);
24 varf vA(u, v) = int3d(cb)(grad(u)' * grad(v)) + on(1, u = 0);
25 varf vB(u, v) = int2d(cb, 2)(u * v);
26 matrix LocA = vA(Vh, Vh);
27 matrix LocB = vB(Vh, Vh);
28 A = LocA;
29 Mat B(A, LocB);
30 int nev = 10;
31 real[int] val(nev);
32 Vh[int] vec(nev);
33 int k = EPSSolve(A, B, vectors = vec, values = val, sparams = "-st_type sinvert -eps_nev " + nev + " -eps_gen_hermitian -eps_view_values ::ascii_matlab");
34 k = min(k, nev);
35 for(int i = 0; i < k; ++i) {
36     macro params()cmm = "Eigenvalue #" + i + " = " + val[i], wait = 1, fill = 1//
37     plotD(cb, vec[i], params);
38     int[int] fforder = [1];
39     savevtk("ev.vtu", cb, vec[i], bin = 1, order = fforder, append = true);
40 }