FreeFEM Documentation on GitHub

stars - forks

Variational Inequality

We present, a classical example of variational inequality.

Let us denote \(\mathcal{C} = \{ u\in H^1_0(\Omega), u \le g \}\)

The problem is:

\[u = arg \min_{u\in \mathcal{C}} J(u) = \frac{1}{2} \int_\Omega \nabla u . \nabla u - \int_\Omega f u\]

where \(f\) and \(g\) are given function.

The solution is a projection on the convex \(\mathcal{C}\) of \(f^\star\) for the scalar product \(((v,w)) = \int_\Omega \nabla v . \nabla w\) of \(H^1_0(\Omega)\) where \(f^\star\) is solution of:

\[(f^\star, v ) = \int_{\Omega}{f v}, \forall v \in H^1_0(`\Omega)\]

The projection on a convex satisfy clearly \(\forall v \in \mathcal{C}, \quad (( u -v , u - \tilde{f} )) \leq 0\), and after expanding, we get the classical inequality:

\[\forall v \in \mathcal{C}, \quad \int_\Omega \nabla(u -v) \nabla u \leq \int_\Omega (u-v) f\]

We can also rewrite the problem as a saddle point problem:

Find \(\lambda, u\) such that:

\[\max_{\lambda\in L^2(\Omega), \lambda\geq 0} \min_{u\in H^1_0(\Omega)} \mathcal{L}(u,\lambda) = \frac{1}{2} \int_\Omega \nabla u . \nabla u - \int_\Omega f u + \int_{\Omega} \lambda (u-g)^+\]

where \(((u-g)^+ = max(0,u-g)\).

This saddle point problem is equivalent to find \(u, \lambda\) such that:

\[\left\{ \begin{array}{cc} \displaystyle \int_\Omega \nabla u . \nabla v + \lambda v^+ \,d\omega= \int_\Omega f u , &\forall v \in H^1_0(\Omega) \cr \displaystyle \int_\Omega \mu (u-g)^+ = 0 , & \forall \mu \in L^2(\Omega) , \mu \geq 0, \lambda \geq 0, \end{array}\right.\]

An algorithm to solve the previous problem is:

  1. k=0, and choose \(\lambda_0\) belong \(H^{-1}(\Omega)\)

  2. Loop on \(k = 0, .....\)

    • set \(\mathcal{I}_{k} = \{ x \in \Omega / \lambda_{k} + c * ( u_{k+1} - g) \leq 0 \}\)

    • \(V_{g,k+1} = \{ v\in H^1_0(\Omega) / v = g\) on \({I}_{k} \}\),

    • \(V_{0,k+1} = \{ v\in H^1_0(\Omega) / v = 0\) on \({I}_{k} \}\),

    • Find \(u_{k+1} \in V_{g,k+1}\) and \(\lambda_{k+1} \in H^{-1}(\Omega)\) such that

      \[\left\{\begin{array}{cc} \displaystyle \int_\Omega \nabla u_{k+1}. \nabla v_{k+1} \,d\omega = \int_\Omega f v_{k+1} , &\forall v_{k+1} \in V_{0,k+1} \cr \displaystyle <\lambda_{k+1},v> = \int_\Omega \nabla u_{k+1}. \nabla v - f v \,d\omega & \end{array}\right.\]

      where \(<,>\) is the duality bracket between \(H^{1}_0(\Omega)\) and \(H^{-1}(\Omega)\), and \(c\) is a penalty constant (large enough).

You can find all the mathematics about this algorithm in [ITO2003] [HINTERMULLER2002].

Now how to do that in FreeFEM? The full example is:

Tip

Variational inequality

 1 load "medit"
 2 
 3 // Parameters
 4 real eps = 1e-5;
 5 real c = 1000; //penalty parameter of the algoritm
 6 real tgv = 1e30; //a huge value for exact penalization
 7 func f = 1; //right hand side function
 8 func fd = 0; //Dirichlet boundary condition function
 9 
10 // Mesh
11 mesh Th = square(20, 20);
12 
13 // Fespace
14 fespace Vh(Th, P1);
15 int n = Vh.ndof; //number of degree of freedom
16 Vh uh, uhp; //u^n+1 and u^n
17 Vh Ik; //to define the set where the containt is reached.
18 Vh g = 0.05; //discret function g
19 Vh lambda = 0;
20 
21 // Problem
22 varf a (uh, vh)
23     = int2d(Th)(
24           dx(uh)*dx(vh)
25         + dy(uh)*dy(vh)
26     )
27     - int2d(Th)(
28           f*vh
29     )
30     + on(1, 2, 3, 4, uh=fd)
31     ;
32 
33 //the mass Matrix construction
34 varf vM (uh, vh) = int2d(Th)(uh*vh);
35 
36 //two versions of the matrix of the problem
37 matrix A = a(Vh, Vh, tgv=tgv, solver=CG); //one changing
38 matrix AA = a(Vh, Vh, solver=CG); //one for computing residual
39 
40 matrix M = vM(Vh, Vh); //to do a fast computing of L^2 norm : sqrt(u'*(w=M*u))
41 
42 real[int] Aiin(n);
43 real[int] Aii = A.diag; //get the diagonal of the matrix
44 real[int] rhs = a(0, Vh, tgv=tgv);
45 
46 // Initialization
47 Ik = 0;
48 uhp = -tgv;
49 
50 // Loop
51 for(int iter = 0; iter < 100; ++iter){
52     // Update
53     real[int] b = rhs; //get a copy of the Right hand side
54     real[int] Ak(n); //the complementary of Ik ( !Ik = (Ik-1))
55     Ak = 1.; Ak -= Ik[];
56     //adding new locking condition on b and on the diagonal if (Ik ==1 )
57     b = Ik[] .* g[]; b *= tgv; b -= Ak .* rhs;
58     Aiin = Ik[] * tgv; Aiin += Ak .* Aii; //set Aii= tgv i in Ik
59     A.diag = Aiin; //set the matrix diagonal
60     set(A, solver=CG); //important to change preconditioning for solving
61 
62     // Solve
63     uh[] = A^-1* b; //solve the problem with more locking condition
64 
65     // Residual
66     lambda[] = AA * uh[]; //compute the residual (fast with matrix)
67     lambda[] += rhs; //remark rhs = -\int f v
68 
69     Ik = (lambda + c*( g- uh)) < 0.; //the new locking value
70 
71     // Plot
72     plot(Ik, wait=true, cmm=" lock set ", value=true, fill=true);
73     plot(uh, wait=true, cmm="uh");
74 
75     // Error
76     //trick to compute L^2 norm of the variation (fast method)
77     real[int] diff(n), Mdiff(n);
78     diff = uh[] - uhp[];
79     Mdiff = M*diff;
80     real err = sqrt(Mdiff'*diff);
81     cout << "|| u_{k=1} - u_{k} ||_2 = " << err << endl;
82 
83     // Stop test
84     if(err < eps) break;
85 
86     // Update
87     uhp[] = uh[];
88 }
89 
90 // Plot
91 medit("uh", Th, uh);

Note

As you can see on this example, some vector, or matrix operator are not implemented so a way is to skip the expression and we use operator +=, -= to merge the result.

Table of content