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: SNESSolve with constraints           *
 8  *******************************************************/
  
10 load "PETSc"
11 macro dimension()2// EOM
12 include "macro_ddm.idp"
 
14 func Pk = P1;
15 func Pkdc = P1dc;
 
17 mesh Th = square(100, 100);
18 Mat A;
19 MatCreate(Th, A, Pk);
 
21 fespace Vh(Th, Pk);
22 fespace Ph(Th, Pkdc);
23 Vh c = 100;
24 real a = 0.001;
25 func real f(real x)   { return (1+a)*x-log(1+x); }
26 func real df(real x)  { return 1+a-1/(1+x); }
27 func real d2f(real x) { return 1/((1+x)*(1+x)); }
 
29 Ph alpha;
30 Ph dfalpha;
31 Ph d2falpha;
32 Vh u;
 
34 varf vRes(w, v) = int2d(Th)(dfalpha*(dx(u)*dx(v) + dy(u)*dy(v))) + on(1, 2, w = u-0);
35 varf vc(w, v) = int2d(Th)(c*v) + on(1, 2, w = 0);
36 varf vJ(w, v) = int2d(Th)(dfalpha*(dx(w)*dx(v) + dy(w)*dy(v)) + d2falpha*(dx(u)*dx(v) + dy(u)*dy(v)) * (dx(u)*dx(w) + dy(u)*dy(w))) + on(1, 2, w = 0);
 
38 func real[int] funcRes(real[int]& inPETSc) {
39     ChangeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
40     alpha = dx(u)*dx(u) + dy(u)*dy(u);
41     dfalpha = df(alpha);
42     d2falpha = 2 * d2f(alpha);
43     real[int] v = vRes(0, Vh, tgv = -1);
44     real[int] outPETSc;
45     ChangeNumbering(A, v, outPETSc);
46     return outPETSc;
47 }
48 func int funcJ(real[int]& inPETSc) {
49     ChangeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
50     matrix Loc = vJ(Vh, Vh, tgv = -1);
51     A = Loc;
52     return 0;
53 }
 
55 real[int] v, b, upper;
56 u[] = vc(0, Vh, tgv = -1);
57 ChangeNumbering(A, u[], b);
58 {
59     Vh upperf = sqrt((x-0.5)^2 + (y-0.5)^2) < 0.4 ? 0.1 : 1.0;
60     plotD(Th, upperf, cmm = "Upper bound");
61     ChangeNumbering(A, upperf[], upper);
62 }
63 v.resize(b.n);
64 v = b;
65 SNESSolve(A, funcJ, funcRes, b, v, xu = upper, sparams = "-snes_view -snes_vi_monitor -snes_type vinewtonrsls -snes_rtol 1.0e-6 -pc_type lu");
66 ChangeNumbering(A, u[], v, inverse = true, exchange = true);
67 plotD(Th, u, cmm = "Solution");