adaptindicatorp2.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

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;
} ;