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: Stokes equation                      *
 8  *******************************************************/
  
10 load "PETSc"
11 macro dimension()2// EOM
12 include "macro_ddm.idp"
 
14 func Pk2 = [P2, P2];
15 func Pk1 = P1;
 
17 mesh Th = square(50, 10, [5 * x, y]);
18 Mat A, C;
19 {
20     DmeshCreate(Th);
21     plotDmesh(Th, wait = 1);
22     {
23         macro def(u)[u, u#B]//
24         macro init(u)[u, u]//
25         MatCreate(Th, A, Pk2);
26     }
27     MatCreate(Th, C, Pk1);
28 }
 
30 macro grad(u)[dx(u), dy(u)]// EOM
31 real Sqrt = sqrt(2.);
32 macro epsilon(u)[dx(u), dy(u#B), (dy(u) + dx(u#B)) / Sqrt]// EOM
33 macro div(u)(dx(u) + dy(u#B))// EOM
 
35 varf vA([u, uB], [v, vB]) = intN(Th)(grad(u)' * grad(v) + grad(uB)' * grad(vB)) + on(1, 3, u = 0, uB = 0) + on(4, u = y*(1.0-y), uB = 0);
36 varf vB([q], [u, uB]) = intN(Th)(-div(u) * q);
37 fespace Vh(Th, Pk2); // local finite element space for the velocities
38 fespace Ph(Th, Pk1); // local finite element space for the pressure
39 matrix LocA = vA(Vh, Vh);
40 matrix LocB = vB(Ph, Vh);
41 Mat B(A, C, LocB);
42 A = LocA;
43 Mat N = [[A , B],
44          [B', 0]];
45 ObjectView(N);
46 set(N, sparams = "-pc_type lu");
47 real[int] rhs(Vh.ndof + Ph.ndof);
48 rhs(0:Vh.ndof - 1) = vA(0, Vh);
49 real[int] sol = N^-1 * rhs;
50 Vh [solX, solY];
51 solX[] = sol(0:Vh.ndof - 1);
52 Ph solP;
53 solP[] = sol(Vh.ndof:sol.n - 1);
54 macro def2(u)[u, u#B]//
55 plotMPI(Th, [solX, solY], Pk2, def2, real, cmm = "Velocities");
56 macro def1(u)u//
57 plotMPI(Th, solP, Pk1, def1, real, cmm = "Pressure");