/usr/share/doc/freefem++/examples/examples++-tutorial/adaptindicatorP2.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 | border ba(t=0,1.0){x=t; y=0; label=1;}; // comment
border bb(t=0,0.5){x=1; y=t; label=2;};
border bc(t=0,0.5){x=1-t; y=0.5;label=3;};
border bd(t=0.5,1){x=0.5; y=t; label=4;};
border be(t=0.5,1){x=1-t; y=1; label=5;};
border bf(t=0.0,1){x=0; y=1-t;label=6;};
mesh Th = buildmesh (ba(6) + bb(4) + bc(4) +bd(4) + be(4) + bf(6));
savemesh(Th,"th.msh");
fespace Vh(Th,P2);
fespace Nh(Th,P0);
Vh u,v;
Nh rho,logrho;
real[int] viso(21);
for (int i=0;i<viso.n;i++)
viso[i]=10.^(+(i-16.)/2.);
real error=0.01;
func f=(x-y);
problem Probem1(u,v,solver=CG,eps=1.0e-6) =
int2d(Th,qforder=5)( u*v*1.0e-10+ dx(u)*dx(v) + dy(u)*dy(v))
+ int2d(Th,qforder=10)( -f*v);
/*
$$\eta_{K} =\left( h_{K}^{2} || f -\Delta u_{{h}} ||_{L^{2}(K)}^{2} +\sum_{e\in \AK} h_{e} \,||\, [ \frac{\partial u_{h}}{\partial n_{k}}] \,||^{2}_{L^{2}(e)} \right)^{\frac{1}{2}}
$$
$$ \rho_{K}= \left( || f -\Delta u_{{h}} ||_{L^{2}(K)}^{2} +\sum_{e\in \AK} \frac{1}{h_{e}} \,||\, [ \frac{\partial u_{h}}{\partial n_{k}}] \,||^{2}_{L^{2}(e)} \right)^{\frac{1}{2}} $$
*/
varf indicator2(uu,chiK) =
intalledges(Th)(chiK*lenEdge*square(jump(N.x*dx(u)+N.y*dy(u))))
+int2d(Th)(chiK*square(hTriangle*(f+dxx(u)+dyy(u))) );
for (int i=0;i< 4;i++)
{
Probem1;
cout << u[].min << " " << u[].max << endl;
plot(u,wait=1);
cout << " indicator2 " << endl;
rho[] = indicator2(0,Nh);
rho=sqrt(rho);
logrho=log10(rho);
cout << "rho = min " << rho[].min << " max=" << rho[].max << endl;
plot(rho,fill=1,wait=1,cmm="indicator density ",ps="rhoP2.eps",value=1,viso=viso,nbiso=viso.n);
plot(logrho,fill=1,wait=1,cmm="log 10 indicator density ",ps="logrhoP2.eps",value=1,nbiso=10);
plot(Th,wait=1,cmm="Mesh ",ps="ThrhoP2.eps");
Th=adaptmesh(Th,[dx(u),dy(u)],err=error,anisomax=1);
plot(Th,wait=1);
u=u;
rho=rho;
error = error/2;
} ;
|