MetricKuate_edp.edp

load "MetricKuate" 
mesh Th=square(5,5,[(x-0.5)*2,(y-0.5)*2]); 
real x0,y0;//  pour definir l'err  forme n lineare en x0,y0 
real coef =1; 
fespace Vh(Th,P1); 
fespace Wh(Th,P2); 
fespace Ph(Th,P0); 
real c=10; 

func f = tanh(c * (sin( (5 * y)) -  (2 * x))) +  (y * x * x) +   pow( y,  3);; 

func fxxx = 0.16e2 * pow(0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1), 0.2e1) * pow(c, 0.3e1) - 0.32e2 * pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1) * (0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1)) * pow(c, 0.3e1); 
func fxxy = -0.40e2 * pow(0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1), 0.2e1) * pow(c, 0.3e1) * cos( (5 * y)) + 0.80e2 * pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1) * (0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1)) * pow(c, 0.3e1) * cos( (5 * y)) + 0.2e1; 
func fxyy = 0.100e3 * pow(0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1), 0.2e1) * pow(c, 0.3e1) * pow(cos( (5 * y)), 0.2e1) - 0.200e3 * pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1) * (0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1)) * pow(c, 0.3e1) * pow(cos( (5 * y)), 0.2e1) - 0.100e3 * tanh(c * (sin( (5 * y)) -  (2 * x))) * (0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1)) * c * c * sin( (5 * y)); 

func fyyy = -0.250e3 * pow(0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1), 0.2e1) * pow(c, 0.3e1) * pow(cos( (5 * y)), 0.3e1) + 0.500e3 * pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1) * (0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1)) * pow(c, 0.3e1) * pow(cos( (5 * y)), 0.3e1) + 0.750e3 * tanh(c * (sin( (5 * y)) -  (2 * x))) * (0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1)) * c * c * cos( (5 * y)) * sin( (5 * y)) - 0.125e3 * (0.1e1 - pow(tanh(c * (sin( (5 * y)) -  (2 * x))), 0.2e1)) * c * cos( (5 * y)) + 0.6e1; 

/* real p=20; real p3=p-3, p321=p*(p-1)*(p-2); func f= x^p + y^p; func fxxx = p321*x; func fyyy = p321*y; func fxxy=0.; func fxyy=0.; */ 
func err=(fxxx*x0*x0*x0+3*fxxy*x0*x0*y0+3*fxyy*x0*y0*y0+fyyy*y0*y0*y0 )*coef; 

for(int i=1;i<4;i++) 
  { 
    Vh m11,m12,m22; 
    coef = 1; 
    Wh f2=f; 
    MetricKuate(Th,200,0.0001,3,err,[m11[],m12[],m22[]],[x0,y0]); 
   // plot(m11,m22,wait=1,cmm="mmmm"); 
    real cc=10; 
    Th=adaptmesh(Th,cc*m11,cc*m12,cc*m22,IsMetric=1,inquire=1,hmin=0.00001,nbvx=1000000); 
    cout << m11[].max << " " << m12[].max << " " << m22[].max << endl; 
    plot(Th,wait=1,ps="Th.eps"); 
    plot(f2,wait=1); 
    Ph eh = (abs(f2-f)); 
    Ph leh= log10(eh); 
    real[int] viso=[-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1]; 
    plot(leh,fill=1,wait=1,viso=viso,value=1,ps="leh.eps"); 
    cout << i << " .... " << eh[].min << " "<< eh[].max << " "<< eh[].sum/eh[].n << " " 
         << int2d(Th)(eh)/Th.area << " " << Th.nt << " " << Th.nv << endl; 
  } 
   Th=square(5,5,[(x-0.5)*2,(y-0.5)*2]); 


 real cerr= 0.005*(0.000961606/0.000582183)^0.66; 
  for(int i=1;i<4;i++) 
  { 
    Vh m11,m12,m22; 
    coef = 1; 
    Wh f2=f; 

    real cc=10; 
    Th=adaptmesh(Th,f,err=cerr,inquire=1,hmin=0.00001,nbvx=1000000); 
    cout << m11[].max << " " << m12[].max << " " << m22[].max << endl; 
    plot(Th,wait=1,ps="Th2.eps"); 
    plot(f2,wait=1); 
    Ph eh = (abs(f2-f)); 
    Ph leh= log10(eh); 
    real[int] viso=[-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1]; 
    plot(leh,fill=1,wait=1,viso=viso,value=1,ps="leh2.eps"); 
    cout << i << " .... " << eh[].min << " "<< eh[].max << " "<< eh[].sum/eh[].n << " " 
         << int2d(Th)(eh)/Th.area << " " << Th.nt << " " << Th.nv << endl; 
  }