load "UMFPACK64" //DefaultSolver=“UMFPACK64”; load "iovtk" verbosity=1; ofstream fout("u.txt"); //ofstream fouti("uex.txt"); ofstream foutou("w.txt"); real K=1.0/40.0; //macro uex(t) (sin(t)*cos(2.*pi*x)*cos(2.*pi*y))// //macro ux(t) (-2.*pi*sin(t)*sin(2.*pi*x)*cos(2.*pi*y))// //macro uy(t) (-2.*pi*sin(t)*cos(2.*pi*x)*sin(2.*pi*y))// real eps=0.2; real w1=3.0; real w2=50.0; real c1=1.0; real c2=1.0; real mu= 1./2.; func real q1 (real t, real s, real w1) { real t1 = t - rint(t); real s1 = s - rint(s); real delta1=1.0/sqrt(w1*pi); if ( square(s1)+square(t1) <= square(delta1)) return w1; else return 0.0; }; func real q2 (real t, real s, real w2) { real t1 = t - rint(t); real s1 = s - rint(s); real delta2=1.0/sqrt(w2*pi); if ( square(s1)+square(t1) <= square(delta2)) return w2; else return 0.0; }; func real q (real t, real s, real w1, real w2){ if ( t < mu) return q1(x/eps,y/eps,w1); else if (t >= mu) return q2(x/eps,y/eps,w2); }; func real qchap (real t, real s,real w2){ if ( t < mu) return 1; else if (t >= mu) return q2(x/eps,y/eps,w2); else return 0.0; }; // initial values for the spatial domain int N=200.; real t=0.0; real dt=0.1; mesh Th=square(N,N,[x,y]); fespace Vh(Th,P1,periodic=[[2,y],[4,y],[1,x],[3,x]]); fespace Wh(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]); Vh uh=0.0; Vh vh; Vh wh=0.0; Vh oldU=0.0; Vh oldW=0.0; real T=1.0;// arbitraire int[int] Order= [1]; string DataName; //Wh uplot; //uplot=q1(x/eps,y/eps,w1); //plot(uplot,fill=true,value=true,ps="q1.eps"); //Wh wplot; //wplot=q2(x/eps,y/eps,w2); //plot(wplot,fill=true,value=true,ps="q2.eps"); Vh fh; func real f (real t) { return K*t*cos(2*pi*x); }; //func real f(real t) //{ //return (cos(t)+8.*pi^2*sin(t))*cos(2.*pi*x)*cos(2.*pi*y)+c1*(-2.*pi*sin(t)*sin(2.*pi*x)*cos(2.*pi*y)) //+c2*(-2.*pi*sin(t)*cos(2.*pi*x)*sin(2.*pi*y))+ q(x,y,w1,w2)*(sin(t)*cos(2.*pi*x)*cos(2.*pi*y)) ; //}; problem het(uh,vh,init=t) = int2d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh))*K*dt) +int2d(Th)((c1*dx(uh)*vh +c2*dy(uh)*vh)*K*dt) +int2d(Th)(K*uh*vh) -int2d(Th)(K*oldU*vh) +int2d(Th)(q(x,y,w1,w2)*uh*vh*dt) -int2d(Th)(fh*vh*dt) ; problem hom(wh,vh,init=t) = int2d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh))*K*dt) +int2d(Th)((c1*dx(wh)*vh +c2*dy(wh)*vh)*K*dt) +int2d(Th)(K*wh*vh) -int2d(Th)(K*oldW*vh) +int2d(Th)(qchap(x,y,w2)*wh*vh*dt) -int2d(Th)(fh*vh*dt) ; for (int n=0;n