/usr/share/doc/freefem++/examples/examples++-tutorial/freeboundary-weak.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 | //
// calcul d'une zone saturation en eau (nappe phréatique)
//
verbosity=1;
real weak=1;
string com= " avec du/dn ";
if(weak) com = " avec du/du = residu ";
real L=10; // longueur du domaine
real Q=0.02; // flux entrant
real K=0.5; //permeabilité
real erradap=0.001;
real coef=1;
real h=2.1; // hauteur du bord gauche
real h1=0.35; // hauteur du bord droite
// maillage d'un tapeze
border a(t=0,L){x=t;y=0;}; // bas
border b(t=0,h1){x=L;y=t;}; // droite
border f(t=L,0){x=t;y=t*(h1-h)/L+h;}; // free surface
border d(t=h,0){x=0;y=t;}; // gauche
int n=10;
mesh Th=buildmesh (a(L*n)+b(h1*n)+f(sqrt(L^2+(h-h1)^2)*n)+d(h*n));
plot(Th,ps="dTh.eps");
fespace Vh(Th,P1);
int j=0,ii=0;
Vh p,v,pp,vv;
varf vPp(p,pp) = int2d(Th)( dx(p)*dx(pp)+dy(p)*dy(pp));
varf vonfree(p,pp) =on(f,p=1);
Vh onfree,wdpdn;
onfree[]=vonfree(0,Vh);
cout << "pb ?" << endl;
onfree[]/=onfree[].max; // sauce
cout << "pb non " << endl;
matrix A= vPp(Vh,Vh,solver=CG);
problem Pp(p,pp,solver=CG) = int2d(Th)( dx(p)*dx(pp)+dy(p)*dy(pp))
+ on(b,f,p=y) ;
// -- imprecise (fort)
problem Pv(v,vv,solver=CG) = int2d(Th)( dx(v)*dx(vv)+dy(v)*dy(vv))
+ on (a, v=0) + int1d(Th,f)(vv*((Q/K)*N.y- (dx(p)*N.x+dy(p)*N.y)));
// -- precise (faible)
problem Pw(v,vv,solver=CG) = int2d(Th)( dx(v)*dx(vv)+dy(v)*dy(vv))
+ on (a, v=0) + int1d(Th,f)(vv*((Q/K)*N.y)) + wdpdn[];
real errv=1;
verbosity=1;
while(errv>1e-6)
{
j++;
Pp;
if (weak) { wdpdn[] = A*p[]; wdpdn[] = wdpdn[].*onfree[]; wdpdn[] = -wdpdn[];
Pw;}
else Pv;
errv=int1d(Th,f)(v*v);
plot(Th,p,v ,cmm=com+" iter = "+j+ " errv =" +errv,wait=0);
// if (j%10==0)
// Th=adaptmesh(Th,p,err=erradap ) ;
real coef=1;
real mintcc = checkmovemesh(Th,[x,y])/5.;
real mint = checkmovemesh(Th,[x,y-v*coef]);
if (mint<mintcc || j%10==0) {
Th=adaptmesh(Th,p,err=erradap ) ;
mintcc = checkmovemesh(Th,[x,y])/5.;
wdpdn=0;
onfree=0;
}
while (1)
{
real mint = checkmovemesh(Th,[x,y-v*coef]);
if (mint>mintcc) break;
cout << " min |T] " << mint << endl;
coef /= 1.5;
}
Th=movemesh(Th,[x,y-coef*v]); // calcul de la deformation
A= vPp(Vh,Vh,solver=CG);
onfree[]=vonfree(0,Vh);onfree[]/=onfree[].max; // sauce
cout << "\n\n"<<j <<"------------ errv = " << errv << "\n\n";
}
//plot(Th,wiat=1,ps="d_Thf.eps",cmm=com);
plot(p,wait=1,ps="d_u.eps",cmm=com);
|