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: Intel MKL PARDISO shared-memory      *
 8  *                parallel solver                      *
 9  *******************************************************/
 
11 load "PARDISO"
12 macro grad(u)[dx(u), dy(u)]// EOM
 
14 func real wedge(real a, real b) {
15     if(y < 0.4 + 0.1 * 0.75 * x)
16         return 2.0;
17     else if(y < 0.8 - 0.2 * 0.75 * x)
18         return 1.5;
19     else
20         return 3.0;
21 }
 
23 mesh Th = square(400, 400);
24 real omega = 2 * pi * 5;
25 func f = 80 * 100 * exp(-20 * 100 * ((x-0.5)^2 + (y-0.25)^2));
26 matrix<complex> A;
27 fespace Ph(Th, P0);
28 Ph val = wedge(x, y);
29 Ph k = omega / val;
30 varf vPb(u, v) = int2d(Th)(-k^2 * u * v + (grad(u)' * grad(v))) + int1d(Th, 2)(1i * k * u * v) + int2d(Th)(f * v) + on(1, u = 0.0);
31 fespace Vh(Th, P2);
32 A = vPb(Vh, Vh);
33 verbosity = 5;
34 set(A, solver = sparsesolver);
35 verbosity = 0;
36 complex[int] b = vPb(0, Vh, tgv = -1);
 
38 Vh<complex> u;
39 u[] = A^-1 * b;
40 plot(u);