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: Mat linear operations                *
 8  *******************************************************/
  
10 load "PETSc"
11 macro dimension()2// EOM
12 include "macro_ddm.idp"
 
14 macro grad(u)[dx(u), dy(u)]// EOM
15 func Pk = P1;
 
17 border upper(t=0, pi)    { x=cos(t); y=sin(t); label=1; }
18 border lower(t=pi, 2*pi) { x=cos(t); y=sin(t); label=2; }
19 mesh Th = buildmesh(upper(100) + lower(100));
20 Mat A;
21 MatCreate(Th, A, Pk);
 
23 fespace Vh(Th, Pk);
24 varf vPb(u, v) = intN(Th)(grad(u)' * grad(v)) + intN(Th)(v) + on(2, u = 0);
25 matrix Loc = vPb(Vh, Vh, tgv = -2);
26 real[int] b = vPb(0, Vh, tgv = -1);
 
28 A = Loc;
 
30 Vh u;
31 set(A, sparams = "-pc_type lu");
32 u[] = A^-1 * b;
33 plotD(Th, u, cmm = "Solution");
 
35 real[int] r = A * u[];
36 exchange(A, b, scaled = true);
37 r -= b;
38 r *= -1;
39 Vh v;
40 v[] = r;
41 plotD(Th, v, cmm = "Residual");
 
43 r = A' * u[];
44 r -= b;
45 r *= -1;
46 v[] = r;
47 plotD(Th, v, cmm = "Residual computed with A'");
 
49 real[int] bPETSc, xPETSc;
50 ChangeNumbering(A, b, bPETSc);
51 xPETSc.resize(bPETSc.n);
52 KSPSolve(A, bPETSc, xPETSc);
53 ChangeNumbering(A, u[], xPETSc, inverse = true, exchange = true);
54 plotD(Th, u, cmm = "Solution with KSPSolve");