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