load “bem”

real k = 10;

int n = 100;

border circle(t = 0, 2*pi){x=cos(t); y=sin(t);}
meshL ThL = buildmeshL(circle(n));
ThL = OrientNormal(ThL,unbounded=1);
varf vbem(u,v) = int1dx1d(ThL)(ThL)(BEM(BemKernel(“SL”,k=k),u,v));
fespace Uh(ThL,P1);
HMatrix H = vbem(Uh,Uh);
func uinc = exp(1i*k

*x);*

Uh p, b;

varf vrhs(u,v) = -int1d(ThL)(uincv);

Uh p, b;

varf vrhs(u,v) = -int1d(ThL)(uinc

b = vrhs(0,Uh);

p = H^-1

*b[];*

varf vpot(u,v) = int1d(ThL)(POT(BemPotential(“SL”,k=k),u,v));

int np = 200;

int R = 4;

border b1(t=-R, R){x=t; y=-R;}

border b2(t=-R, R){x=R; y=t;}

border b3(t=-R, R){x=-t; y=R;}

border b4(t=-R, R){x=-R; y=-t;}

mesh ThOut = buildmesh(b1(np)+b2(np)+b3(np)+b4(np)+circle(-n));

fespace UhOut(ThOut,P1);

HMatrix HP = vpot(Uh,UhOut);

UhOut u, utot;

u[] = HPp;

varf vpot(u,v) = int1d(ThL)(POT(BemPotential(“SL”,k=k),u,v));

int np = 200;

int R = 4;

border b1(t=-R, R){x=t; y=-R;}

border b2(t=-R, R){x=R; y=t;}

border b3(t=-R, R){x=-t; y=R;}

border b4(t=-R, R){x=-R; y=-t;}

mesh ThOut = buildmesh(b1(np)+b2(np)+b3(np)+b4(np)+circle(-n));

fespace UhOut(ThOut,P1);

HMatrix HP = vpot(Uh,UhOut);

UhOut u, utot;

u[] = HP

utot = u + uinc;

plot(utot,fill=1,value=1,cmm=“u_total”);