include "getARGV.idp" real x0=-7.,y0=-7.; real xF=7.,yF=7.; real[int] colortab(-1.1:0.1:1.1); int n=7; int N =2^n; real h=(xF-x0)/N; real xS=-1.4,yS=-1.4; real xSF=1.4,ySF=1.4; 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,T1; dt=(0.163*h);// step in time T1=20.; td=floor(T1/dt); real idt2=1/(dt^2); real idt1=1/(2*dt); real theta=pi/2.; real lambdahalf=1.4; real xbeam=xF+lambdahalf*cos(theta),ybeam=yF+lambdahalf*sin(theta); // domain of problem int[int] outer=[1,1,1,1]; mesh Th = square(N,N,[x0+(xF-x0)*x ,y0+(yF-y0)*y],label=outer); int[int]labs1=[0,1]; Th=change(Th,region=labs1);//chang the region number from 0 to 1 int[int] inner=[2,2,2,2]; mesh ThS = square(N,N,[xS+(xSF-xS)*x ,yS+(ySF-yS)*y],label=inner);//scatterd region int[int]labs2=[0,2]; ThS=change(ThS,region=labs2);// change thr region number from 0 to 2 int outerboundary=1; mesh Thwithole=trunc(Th,(xxSF|yxSF| y>ySF ),label=outerboundary); Thwithole=change(Thwithole,region=labs1); int innerboundary=2; Thwithole=trunc(Th,(xxSF|yxSF| y>ySF ),label=innerboundary); int Rlabel1= Th(7.0,7.0).region; int Rlabel2=ThS(1.4,1.4).region; int Rlabel3=Thwithole(7.0,7.0).region; int Rlabel4=Thwithole(1.4,1.4).region; /* cout<<"Region label in point (7.0, 7.0) in Th =" << Rlabel1<-1 && zz<1))//EOM macro utilde(zz) ((-1)*h1(z0(x,y)))//EOM problem wave(u,v)=int2d( Thwithole)(idt2*u*v)-int2d( Thwithole)(2*idt2*uold*v)+int2d ( Thwithole)(idt2*uvold*v)+int2d( Thwithole)(c^2*grad(uold)'*grad(v)) +int1d(Thwithole,outer)((c*u*v)*idt1)-int1d( Thwithole,outer)((c*uvold*v)*idt1) +on(inner, u=utilde(zz)); // loop on the time for a given mesh for(int s=0;s<=td;s++) { T =s*dt; //initial conditions verbosity = getARGV("-verbosity", 0); wave; uvold=uold; uold =u; plot(Thwithole,u,cmm= "Forward problem Wave equation u_Scattered field domain with hole t="+T+", N = "+N,fill=1,viso=colortab,dim=2,value=true); ofstream s("u_S with hole="+s+" for N="+N+".txt"); s<