/usr/share/doc/freefem++/examples/examples++/Richard.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 126 | real Ks=0.01,
hg=30,
thetas=0.3,
eta = 6.55,
m = 0.173,
n = 2/(1-m);
real z0=215;
real q0=15/3600.;
real dt=60;
real hmax=0;
// $A(h) - \partial h / \partial t - div(K(h)(\nabla(h-y)) = f $ dans $ \Omega$
// -K(h)(\nabla(h-y)). n = q0 $ sur $ \Gamma_1$
// h = h_0$ sur $\Gamma_0$
// + condition initial $d_d$
// A(h) = h <0 ? C(h) : 0;
//
// remarque z == y
real xmax = 300, ymax=300, x0=60, y0= 215;
border ba(t=0,ymax) { x=0; y=ymax-t ;label=2;}; // gauche
border bb1(t=x0,0) { x=t; y=ymax ;label=1;}; // haut 1
border bb2(t=xmax,x0) { x=t; y=ymax ;label=2;}; // haut 2
border bc1(t=y0,0) { x=xmax; y=ymax-t ;label=2;}; // droite
border bc2(t=ymax,y0) { x=xmax; y=ymax-t ;label=3;}; // droite
border bd(t=0,xmax) { x=t; y=0; label=4;}; // bas
int Gamma0=3;
int Gamma1=1;
int nn=20;
int nn1=nn*x0/xmax,nn2=nn-nn1;
int ny1=nn*y0/ymax,ny2=nn-ny1;
plot(ba(nn)+bb1(nn1)+bb2(nn2)+bc1(ny1)+bc2(ny2)+bd(nn),wait=1);
mesh Th=buildmesh(ba(nn)+bb1(nn1)+bb2(nn2)+bc1(ny1)+bc2(ny2)+bd(nn));
plot(Th,wait=1);
fespace Vh(Th,P1);
Vh h,v,hhh;
macro theta(h) (thetas*(1+((abs(h)/hg))^n)^(-m))//
macro dtheta(h) (m*n*thetas*(1+((abs(h)/hg))^n)^(-m-1)*(((abs(h)/hg))^(n-1))/hg)
//
macro A(h) ( (h<=0)* dtheta(h) )
//
macro K(h) (Ks*((h<=0)*((theta(h)/thetas)^eta)+ (h>0)))
//
Vh hd= -y0+(ymax-y); // bof bof ????
Vh hn=hd,hh;
Vh Ahdt,Kh;
int nbiso=20;
real[int] viso(3+(75+110/2)/5);
{int k=0;
for(int i=-75;i<0;i+=5)
viso(k++)=i;
viso(k++)=-0.5;
viso(k++)=0.;
viso(k++)=0.5;
for(int i=5;i<=110;i+=5*2)
viso(k++)=i;
}
/*
problem Richard(h,v) =
int2d(Th)( Ahdt * h * v+ Kh* (dx(h)*dx(v)+dy(h)*dy(v)) )
- int2d(Th)( Ahdt* hn*v - Kh* dy(v) )
- int1d(Th,Gamma1)(q0*v)
+ on(3,h=(ymax-y)-y0)
;
*/
real pena=1e10;
problem Richard(h,v) =
int2d(Th)( Ahdt * h * v+ Kh* (dx(h)*dx(v)+dy(h)*dy(v)) )
- int2d(Th)( Ahdt* hn*v - Kh* dy(v) )
- int1d(Th,Gamma1)(q0*v)
+int1d(Th,3)(pena*h*v)-int1d(Th,3)(pena*((ymax-y)-y0)* v)
;
plot(hn,wait=1,cmm=" hd ");
Ahdt=0; Kh=1;
// Richard;
// plot(hd,wait=1,cmm="hd ----");
real temps=0;
for(int i=0;i<100;i++)
{
string scmm="h + temps "+int(temps)/3600+"h "+ ((temps)%3600/60.) + "mn ";
for(int k=0;k<3;k++)
{
Kh=K(h);
Ahdt=A(h)/dt;
cout << " "<< Kh[].min << " " << Kh[].max << endl;
// plot(Ahdt,fill=1,value=1,wait=1,cmm="Ahdt");
// plot(Kh,fill=1,value=1,wait=1,cmm="Kh");
// plot(Kh,wait=1,cmm="Kh");
Richard;
cout << " h: min " << h[].min << " max " << h[].max <<endl;
hmax=h[].max ;
// plot(h,wait=0,cmm=scmm,viso=viso);
// hhh = h <0;
}
if(i%10==1) {
Th=adaptmesh(Th,h,ratio=1.1);
// plot(Th,h,cmm="h ",value=1,wait=1);
}
// plot(hhh,wait=1,cmm="h < 0");
plot(h,wait=0,cmm=scmm,viso=viso,value=1);
// plot(h,cmm="h ",value=1);
hn=h;
temps += dt;
}
cout << " hmax = " << hmax << endl;
|