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: Liouville–Bratu–Gelfand equation     *
 8  *******************************************************/
  
10 load "PETSc"
11 macro dimension()2// EOM
12 include "macro_ddm.idp"
 
14 real c = 6.25;
15 int N  = 80;
 
17 mesh M;
18 {
19     border a(t=0,1) { x=t;   y=0;   label=1; };
20     border b(t=0,1) { x=1;   y=t;   label=1; };
21     border c(t=0,1) { x=1-t; y=1;   label=1; };
22     border d(t=0,1) { x=0;   y=1-t; label=1; };
23     M = buildmesh(a(N)+b(N)+c(N)+d(N));
24 }
 
26 load "Element_P3"
27 func Pk = P3;
28 Mat J;
29 MatCreate(M, J, P3);
30 fespace Vh(M, Pk);
31 Vh u;
32 func BC = cos(pi*x)*cos(pi*y);
33 varf vInit(w, v) = on(1, w = BC);
34 varf vJ(w, v) = int2d(M)(dx(w)*dx(v) + dy(w)*dy(v) - c*exp(u)*w*v) + on(1, w = 0);
35 varf vRes(w, v) = int2d(M)(dx(u)*dx(v) + dy(u)*dy(v) - c*exp(u)*v) + on(1, w = u);
36 int it = 20;
37 real tol = 1.0e-6;
38 u[] = vInit(0, Vh, tgv = -1);
39 real[int] b = u[];
40 for(int i = 1; i < it; ++i) {
41     real[int] v = vRes(0, Vh, tgv = -1);
42     v -= b;
43     real res = sqrt(J(v, v));
44     if(mpirank == 0)
45         cout << "Norm of the residual = " << res << endl;
46     if(res < tol)
47         break;
48     matrix Loc = vJ(Vh, Vh, tgv = -1);
49     J = Loc;
50     real[int] inc = J^-1 * v;
51     u[] -= inc;
52     plotD(M, u, cmm = "Solution at iteration #" + i);
53 }
54 plotD(M, u, cmm = "Final solution");