Actual source code: example2.edp
 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: 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 createMat(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")