1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
// variationnal inequality
// --------------------------
// Probleme:
// $ - \Delta u = f , \quad u=g \on \Gamma, \quad u < umax $
// algo of Primal-Dual Active set strategy as a semi smoth Newton Method
// HinterMuller , K. Ito, K. Kunisch
// to appeared in SIAM Option
// Thank to O. Pironneau
// --------------------------
// F. Hecht
// -----------------------
mesh Th=square(10,10);
real eps=1e-10;
fespace Vh(Th,P1); // P1 FE space
int n = Vh.ndof; // number of Degree of freedom
Vh uh,uhp,Ah; // unkown and test function.
real[int] rhs(n);
real cc=1000;
func f=1; // right hand side function
func g=0; // boundary condition function
func gmax=0.05;
Vh umax=gmax;
//
real tol=0.05,tolmin=0.002;
real tgv = 1e30;
real res=0;
varf a(uh,vh) = // definion of the problem
int2d(Th)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) // bilinear form
- int2d(Th)( f*vh ) // linear form
+ on(1,2,3,4,uh=g) ; // boundary condition form
varf vM(uh,vh) = int2d(Th)(uh*vh);
matrix A=a(Vh,Vh,tgv=tgv,solver=CG);
matrix AA=a(Vh,Vh);
matrix M=vM(Vh,Vh);
rhs = a(0,Vh,tgv=tgv);
real[int] Aii(n),Aiin(n),Ah1(n),b(n);
Aii=A.diag; // get the diagonal of the matrix
//cout << " Aii= " << Aii << endl;
Ah =0;
uhp=0;
Vh lh=0;
int kadapt=0,kkadapt=0;
for(int iter=0;iter<100;++iter)
{
// solve the problem plot(uh); // to see the result
b = rhs;
// add new lock condition on i / if (Ah[i] ==1 )
Ah1= 1.; Ah1 -= Ah[]; // Ah1 = ! Ah
b = Ah[] .* umax[]; b *= tgv; b -= Ah1 .* rhs;
Aiin = Ah[] * tgv; Aiin += Ah1 .* Aii;
A.diag = Aiin;
set(A,solver=CG); // important to change precondiconning
uh[] = A^-1* b;
lh[] = AA * uh[];
lh[] += rhs;
// plot(lh,wait=1);
Ah = ( lh + cc*( umax- uh)) < 0.;
// plot(Ah, wait=1,cmm=" lock ",value=1 );
plot(uh,wait=1,cmm="uh");
real[int] d(n),Md(n);
d= uh[]-uhp[];
Md = M*d;
real err = sqrt(Md'*d);
Md=M*uh[];
Ah1=1.;
real intuh = (Ah1'*Md); // int uh;
cout << " err norm L2 " << err << " "
<< " int uh = " << intuh
<< " kkadapt =" << kkadapt <<endl;
res = intuh;
if(err< eps && kkadapt ) break;
bool adapt = err< eps || (iter%5 == 4) ;
if(adapt)
{
kadapt++;
Th=adaptmesh(Th,uh,err=tol);
kkadapt = tol == tolmin; // we reacht the bound
tol = max(tol/2,tolmin);
cout << " ++ tol = " << tol << " " << kadapt << " " << kkadapt <<endl;
plot(Th,wait=0);
uhp=uhp;
n=uhp.n;
uh=uh;
Ah=Ah;
lh=lh;
umax = gmax;
A=a(Vh,Vh,tgv=tgv,solver=CG);
AA=a(Vh,Vh);
M=vM(Vh,Vh);
Aii.resize(n);
Aiin.resize(n);
Ah1.resize(n);
b.resize(n);
rhs.resize(n);
Aii=A.diag; // get the diagonal of the matrix
rhs = a(0,Vh,tgv=tgv);
}
uhp[]=uh[];
}
savemesh(Th,"mm",[x,y,uh*10]);
{
int nn=100;
real[int] xx(nn+1),yy(nn+1),pp(nn+1);
for (int i=0;i<=nn;i++)
{
xx[i]=i/real(nn);
yy[i]=uh(0.5,i/real(nn));
}
plot([xx,yy],wait=1,cmm="u1 x=0.5 cup");
}
Aiin=M*uh[];
Aii=1.;
cout << " -- int uh = " << res << endl;
assert( abs(res-0.0288611) < 0.001);
|
