include "getARGV.idp" real x0=0.0,y0=0.0; real xF=3.0,yF=3.0; int n=7; int N =2^n; real h=(xF-x0)/N; real xI=0.15,yI=0.15; real xIF=0.75,yIF=0.75; real xmid=(xI+xIF)*0.5; real ymid=(yI+yIF)*0.5; real xS=1.05,yS=1.05; real xSF=1.95,ySF=1.95; func c1=1e-4*((x>xS)*(xyS)*(yxSF|yxSF| y>ySF ));// velocity of the wave without obstacle func c=c1+c2; macro grad(u) [dx(u), dy(u)]//EOM real t,dt,td,T; dt=0.163*h;// step in time td=floor(3/dt); real idt2=1/(dt^2); real idt1=1/(2*dt); int[int] Energy(786); // domain of problem mesh Th = square(N,N,[x0+(xF-x0)*x,y0+(yF-y0)*y]); int[int] innerboundary=[1,1,1,1]; mesh ThS=square(N,N,[xS+(xSF-xS)*x,yS+(ySF-yS)*y],label=innerboundary );// subdomain for the obstacle fespace Vh(Th,P1); Vh u,v,uold,uvold,v0,g,g1,k,k1,H,H1,f; fespace Ph(Th,P0); Ph cvelocity = c; plot(cvelocity,cmm=" velocity in the domain", wait=1, fill=1,value=1); g=(1+2.*(xmid-x)/(xmid-xI))*((xI-x)/(xI-xmid))^2*(x>=xI)*(x<=xmid); k=(1+2.*(x-xmid)/(xIF-xmid))*((xIF-x)/(xIF-xmid))^2*(x>=xmid)*(x<=xIF); H=g+k; g1=(1+2.*(ymid-y)/(ymid-yI))*((yI-y)/(yI-ymid))^2*(y>=yI)*(y<=ymid); k1=(1+2.*(y-ymid)/(yIF-ymid))*((yIF-y)/(yIF-ymid))^2*(y>=ymid)*(y<=yIF); H1=g1+k1; // initialisation f=H*H1;//Hermit cubic bell on ThI uvold=f; uold=uvold; plot(Th,uvold,cmm= "initial time u0",wait=true,fill=1,dim=2,value=true); // loop on the time for a given mesh for(int s=0;s<=td;s++) { T =s*dt; problem wave(u,v)=int2d(Th)(idt2*u*v)-int2d(Th)(2*idt2*uold*v)+int2d(Th)(idt2*uvold*v) +int2d(Th)(c^2*grad(uold)'*grad(v)) +on(1,2,3,4,u=0); //initial conditions verbosity = getARGV("-verbosity", 0); wave; Energy[s]=int2d(Th)(((u-uold)/dt)^2+c^2*grad(u)'*grad(u)); cout<<""<