include "getARGV.idp" real x0=0.0,y0=0.0; real xF=3.0,yF=3.0; int N =2^7; real h=(xF-x0)/N;//step size in mesh real xS=1.05,yS=1.05; real xSF=1.95,ySF=1.95; real xTRAC=0.9,yTRAC=0.9; real xTRACF=2.1,yTRACF=2.1; func c1=1e-4*((x>xS)*(xyS)*(yxSF|yxSF| y>ySF ));// velocity of the wave without the obstacle //define the sensors points real x1=3,y1=8*h; real x2=3,y2=24*h; real x3=3,y3=40*h; real x4=3,y4=56*h; real x5=3,y5=72*h; real x6=3,y6=88*h; real x7=3,y7=104*h; real x8=3,y8=120*h; macro grad(u) [dx(u), dy(u)]//EOM real t, dt,td,T; dt=0.163*h;// step in time real idt2=1/(dt^2); real idt1=1/(2*dt); real r=h/4; func c=c1+c2; td=floor(3/dt); // domain of problem int[int] lab = [1,2,10,5]; mesh ThBot = square(N,N,[x0+(xF-x0)*x,y0+0.5*(yF-y0)*y] , label = lab); lab = [10,2,3,4]; mesh ThTop = square(N,N,[x0+(xF-x0)*x,y0+0.5*(yF-y0)*y+1.5], label = lab); mesh Th = ThBot + ThTop; Th = change(Th, rmInternalEdges = true); int[int] bTRAC=[1,1,1,1]; mesh ThTRAC = square(N,N,[xTRAC+(xTRACF-xTRAC)*x,yTRAC+(yTRACF-yTRAC)*y], label = bTRAC); mesh Thm=Th+ThTRAC; Thm = change(Thm, rmInternalEdges = true); int outboundary=5; mesh Thwithole= trunc(Thm,((xxSF|yxSF| y>ySF )),label= outboundary); int innerboundary=2; Thwithole= trunc(Thm,((xxSF|yxSF| y>ySF )),label= innerboundary); plot(Thwithole,wait=1); /* fespace Vh(Thwithole,P1); Vh u,v,uold,uvold,uoldtemp,kappa1,kappa2,kappa3,kappa4,kappa5,kappa6,kappa7,kappa8,r1,p; //define the sensors region- eight sensors only border b1(t1=0,2*pi){x=x1+r*cos(t1); y=y1+r*sin(t1);} mesh Th1=buildmesh(b1(3));//three nodes mesh eTh1 = emptymesh(Th1); kappa1=(chi(eTh1)(x,y))*((x-x1)^2+(y-y1)^2<=r^2); border b2(t1=0,2*pi){x=x2+r*cos(t1); y=y2+r*sin(t1);} mesh Th2=buildmesh(b2(3));//three nodes mesh eTh2 = emptymesh(Th2); kappa2=(chi(eTh2)(x,y))*((x-x2)^2+(y-y2)^2<=r^2); border b3(t1=0,2*pi){x=x3+r*cos(t1); y=y3+r*sin(t1);} mesh Th3=buildmesh(b3(3));//three nodes mesh eTh3 = emptymesh(Th3); kappa3=(chi(eTh3)(x,y))*((x-x3)^2+(y-y3)^2<=r^2); border b4(t1=0,2*pi){x=x4+r*cos(t1); y=y4+r*sin(t1);} mesh Th4=buildmesh(b4(3));//three nodes mesh eTh4 = emptymesh(Th4); kappa4=(chi(eTh4)(x,y))*((x-x4)^2+(y-y4)^2<=r^2); border b5(t1=0,2*pi){x=x5+r*cos(t1); y=y5+r*sin(t1);} mesh Th5=buildmesh(b5(3));//three nodes mesh eTh5 = emptymesh(Th5); kappa5=(chi(eTh5)(x,y))*((x-x5)^2+(y-y5)^2<=r^2); border b6(t1=0,2*pi){x=x6+r*cos(t1); y=y6+r*sin(t1);} mesh Th6=buildmesh(b6(3));//three nodes mesh eTh6 = emptymesh(Th6); kappa6=(chi(eTh6)(x,y))*((x-x6)^2+(y-y6)^2<=r^2); border b7(t1=0,2*pi){x=x7+r*cos(t1); y=y7+r*sin(t1);} mesh Th7=buildmesh(b7(3));//three nodes mesh eTh7 = emptymesh(Th7); kappa7=(chi(eTh7)(x,y))*((x-x7)^2+(y-y7)^2<=r^2); border b8(t1=0,2*pi){x=x8+r*cos(t1); y=y8+r*sin(t1);} mesh Th8=buildmesh(b8(3));//three nodes mesh eTh8 = emptymesh(Th8); kappa8=(chi(eTh8)(x,y))*((x-x8)^2+(y-y8)^2<=r^2); //weak formulation problem backwardwave(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, bTRAC)((c*u*v)*idt1)-int1d(Thwithole,bTRAC)((c*uvold*v)*idt1) +on (outboundary,u=0); verbosity = getARGV("-verbosity", 0); ifstream s785("Results wave at iteration=785 for N=128.txt"); s785>>r1[]; uvold=r1*(kappa1+kappa2+kappa3+kappa4+kappa5+kappa6+kappa7+kappa8);//u_T in the sensors ifstream s784("Results wave at iteration=784 for N=128.txt"); s784>>p[]; uold=p*(kappa1+kappa2+kappa3+kappa4+kappa5+kappa6+kappa7+kappa8);//u_T-1 in the sensors for(int j=2;j<=td;j++) { int k=td-j; T=k*dt; Vh rk; ifstream sk("Results wave at iteration="+k+" for N="+N+".txt"); sk>>rk[]; backwardwave; uvold=uold; uold=u; uold=uold+rk*(kappa1+kappa2+kappa3+kappa4+kappa5+kappa6+kappa7+kappa8); plot(Thwithole,u,cmm="backward problem TRAC T="+T+"",fill=1,dim=2,value=true); } */