/usr/share/doc/freefem++/examples/examples++-tutorial/VI-adap.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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);
|